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ABSTRACT 


A general, modular, pulse propagation model for underwater acoustics that is based on 
linear systems theory for sound-speed profiles as a function of depth is presented. The 
development and computer implementation of the model, together with results from 
preliminary computer simulation studies involving the transmission of CW and LFM pulses 
from a planar array of complex weighted point sources is reported. The studies examined 
free-space propagation problems (i.e., no boundaries) in homogeneous and 
inhomogeneous media using a transfer function of the ocean medium based on the WKB 
approximation. The two main outputs from the model are the predicted complex acoustic 
field as a function of frequency and spatial location and the time-domain output electrical 


signal from each element in a receive planar array. 
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I. INTRODUCTION 


For small-amplitude acoustic signals it can be shown that the linear acoustic wave 
equation governs the propagation of sound in a fluid medium [Ref. 1:pp. 23-31] (Ref. 
2:pp. 98-105]. Due to the linear nature of this equation, the ocean medium can be treated 
in general, as a linear, time-variant, space-variant, random filter or communication channel 
completely characterized by a transfer function [Ref. 3:pp. 1-6]. As a result, linear 
systems theory can be applied to the ocean medium propagation problem and through the 
use of coupling equations (Ref. 4] analytical expressions can be derived for the complex 
acoustic field and the output electrical signal at each element in a hydrophone array in terms 


of : 


¢ the frequency spectrum of the transmitted pulse 
¢ the far-field directivity function or beam pattern of the transmit array 
¢ the ocean medium transfer function 


¢ the far-field directivity function of the receive array. 


These expressions are developed for planar arrays, since they include linear arrays and 
single transducers as special cases, and are amenable to computer simulation studies. 

The application of Fourier analysis and linear systems concepts to problems in physics 
is not unique and has been applied extensively to the study of optics [Ref. 5]. What is 
unique is the application of these ideas to the physics of ocean acoustics and the treatment 
of the ocean medium, in general, as a time-variant, space-variant linear system [Ref. 3:pp. 


7-28] (Ref. 6]. 


In this thesis, we will present a general, modular, pulse propagation model for 
underwater acoustics that is based on linear systems theory and the coupling equations for 
sound-speed profiles that are a function of depth. Since the coupling equations are the 
formal solution of the pulse propagation problem, and since the coupling equations depend 
on the transfer function of the ocean medium, the only thing that changes from problem to 
problem from an ocean acoustics point of view, is the ocean medium transfer function. 
Therefore, regardless of the problem under consideration, the coupling equations only need 
to be programmed once. Preliminary computer simulation studies can then carried out for 
the specific case of a free-space propagation problem using a transfer function based on the 
WKB approximation. 

The evaluation of the ocean medium transfer function is based on the full-wave solution 
technique, which involves the evaluation of wave propagation vector component integrals. 
This method is different from, but can be related to, ray acoustics and normal mode 
approaches to solving the propagation problem. As a consequence, our model of the pulse 
propagation problem is termed a full-wave solution. The computer implementation of a 
full-wave solution to the propagation problem is not new [Ref. 7] but an implementation of 
a model also based on linear systems concepts is. 

For the purposes of this thesis the transmit and receive arrays were assumed to be 
motionless which meant that we were able to treat the ocean medium as a linear, time- 
invariant, space-variant filter. Chapter 2 introduces and applies the linear systems 
approach to the time-invariant ocean medium problem and develops the coupling equations 
as far as possible without making any further simplifying assumptions about the ocean 
medium. As a consequence the results obtained can be applied to the propagation of an 
acoustic signal in any fluid medium. Chapter 3 then develops the coupling equations for 


the case of an ocean medium that is axially symmetric about the depth axis. This 


development was based on the fact that axial symmetry was due to the speed of sound 
being a function of depth only, i.e., the case we wish to solve for. The results of this 
analysis can then be applied to any ocean medium description that exhibits axial symmetry. 
Chapter 4 then goes on to develop the ocean medium transfer function for the specific case 
of the ocean modelled with the WKB approximation. 

At this point, the development of the pulse propagation model is complete and is 
capable of producing two major outputs. The first output provides the magnitude and 
phase of the complex acoustic field incident upon the receive array as a function of 
frequency and spatial coordinates. This information is important, for example, for target 
localization using matched-field techniques. The second output is the time-domain output 
electrical signal at each element in the receive array. This information is important, for 
example, to illustrate pulse distortion due to the effects of dispersion in a waveguide, and 
as input to space-time signal processing algorithms. 

An integral part of the model was the development, in Chapter 5, of a signal generator 
to simulate the transmitted electrical pulses that would be converted to acoustic energy for 
the propagation problem. It was assumed that the pulses were narrowband amplitude-and 
angle-modulated carriers. The signal generation scheme was based on a truncated Fourier 
series and complex envelope representation of narrowband signals. 

Chapters 6 and 7 contain computer simulation results for the specific case of the 
ocean medium that is characterized by the WKB approximation. These chapters deal with 
simulation results for an ocean medium characterized by a constant speed of sound and a 


sound-speed profile with a single constant gradient, respectively. 


Il. LINEAR SYSTEMS 


A. THE LINEAR SYSTEMS APPROACH 

Treating acoustic signals as small fluctuating disturbances in a fluid medium, it can be 
shown that the principles of fluid mechanics can be used to derive an equation governing 
the propagation of sound in a fluid medium. The linear acoustic wave equation, given by 


2 be ae 
V @(ur) ->—— O(br) = xm (tn), (2.1) 
c (r) ot 


(where @(t,r) is the velocity potential, in m2 / sec, at time t and spatial location r, xj (t,r) 
is the input acoustic signal to the medium, in 1 / sec, and c(r) is the speed of sound in the 
medium, in m/ sec) can be derived by making a number of restrictive assumptions and 
applying them to the equations of state, continuity and momentum for a fluid. These 
restrictions allow us to derive a simple equation, of the form given by Eq. (2.1), which 
nonetheless is adequate to describe most commonly encountered acoustic phenomena. 
[Ref. 1:pp. 23-31] [Ref. 2:pp. 98-105] [Ref. 3:pp. 1-3] 

The more commonly encountered acoustic quantities of acoustic pressure, p(t,r), in 
N / m2 and the acoustic fluid particle velocity, u(t,r), in m / sec can be obtained from the 


velocity potential as follows : 


p(tr) =-p, 2 p(t,r) (2.2) 


where Pq is the constant equilibrium density of the fluid in kg / m3, and 


u(t,r) = Vo(tr). (2.3) 


It should be noted here that the concept of a scalar velocity potential arises from one of the 
restrictions used in the development of Eq. (2.1), namely that the particle motions 
associated with sound waves are irrotational. Therefore, the relationships given by Eqs. 
(2.2) and (2.3) are governed by this assumption. 

The linear acoustic wave equation given by Eq. (2.1) is so called because it is a linear 
partial differential equation relating the velocity potential and input signal. Since Eq. (2.1) 


can be thought of as a linear operation denoted by 


Re 


L{s} = V a a 7 
c(r) ar 





(2.4) 


that relates the velocity potential and input signal, we can use a linear systems theory 
approach to describe the effect of the fluid medium on the propagation of an input signal 
through the medium. 

Using a linear systems approach, we describe the fluid medium as being a linear system 
that is characterized by an impulse response that is a function of both time and space. To 
then find the output (i.e., the acoustic velocity potential @(t,r)) from such a system due to 
some input (i.e., the acoustic signal x, (t,r)), we simply convolve the input with the 
impulse response of the system. The basic geometry of this problem is illustrated in 
Ere, 2.1. 

This operation can also be done in the frequency (in Hz ) and spatial frequency (in 


cycles /m) domains using the transfer function of the system. Since this is analogous to 
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Fig. 2.1 Ocean Medium Geometry 


the concept of filtering, we say that our system is being modelled as a time-variant, space- 
variant, linear filter. [Ref. 3:pp. 3-6] 

Up until this point, we have been talking about any general fluid medium. We now 
wish to consider this problem with specific reference to the ocean medium. For the 
purposes of this thesis, we will model the ocean medium as a time-invariant, space-variant 
linear filter. Time invariance means that we will not be considering Doppler shifts between 
the input and output of the filter. In relation to the ocean medium we are modelling, this 
means that we are ignoring such things as motion between the transmitter and receiver, 
target motion and discrete point scatterers. 

Space variance means that there is a shift in spatial frequency between the transmitted 
and received signal. This means that sound speed is a function of position and there exists 


a non-uniform sound-speed profile. As a result, a transmit signal will scatter due to 


refraction. If we think of a ray interpretation of acoustic propagation, the rays will bend 


rather than travel in straight lines as is the case for homogeneous media. 


B. THE LINEAR, TIME-INVARIANT, SPACE-VARIANT FILTER 


Consider a system modelled by a linear, time-invariant, space-variant filter. The time- 


invariant, space-variant impulse response of the filter, h(t,rg;r), will then completely 


characterize the system being modeled. The relationship between the input x(t,r) and the 


output y(t,r) of such a filter is given by 
y(t,r) -{ f X(t-T,F-Ig )h(t,rg;r) dt dry (2.5) 


where the function h(t,r;r) is the response of the filter positioned at spatial location 
r = (x,y,z), due to an impulse applied at position fg = (X9,Yo.Z9), tT Seconds ago. For 
the specific case of modeling the ocean medium, the terms y(t,r) and x(t,r) in Eq. (2.5) 
would equate to 9(t,r) and Xyyq (t.r) in Fig. 2.1, respectively. 
As an input to our filter, consider a signal x(t,r) of the following form: 

Raise ee (2.6) 
where f is the input frequency in Hz and v = (fy ,fy,fz) are the input spatial frequency 
components in cycles / meter. The spatial frequencies are a function of direction and 


wavelength given by 


ah (2.7a) 


fy =v/A (2.7b) 


and 


fy =w/d (2.7c) 


where u, v and w are the direction cosines with respect to the positive X, Y and Z axes. 
The surfaces of constant phase, commonly referred to as wavefronts, for the signal 
described by Eq. (2.6) are given by 2aver = constant, i.e. a plane., Consequently, Eq. 
(2.6) describes a time-harmonic plane wave traveling in the direction of the propagation 
vector k = 21v, normal to the wavefront.[Ref. 2: pp. 107-108] 

To find the filter output due to a time-harmonic plane wave input, we combine 
Eqs. (2.5) and (2.6) which yields 


j2nft -j2mve 
jonft -j2nver 


y(t,r) =e H(f,v;r), (2.8) 


where 


H(f,vsr) = f f h(tyroir) eo eo dt dr, (2.9) 


is the transfer function of the linear, time-invariant, space-variant system (or filter). It is 
completely analogous to the transfer function H(f) of linear time-invariant (LTI) systems 


often encountered in electrical circuit theory and signal analysis. 


Equation (2.9) is in the form of a multidimensional Fourier transform of the impulse 


response and can be written as 


HG) ibs F, { h(tnose (2.10) 
where F, represents the forward time-domain Fourier transform given by 


H(f,ro:r) =F, ( h(trgsr) ) = f h(tyrow)e dt, (2.11) 
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and ry represents the forward spatial Fourier transform given by 


Hvir) =F, (h(tro;r)} = f h(trgir) Yo dre . (2.12) 


Oo 


Comparing Eqs. (2.8) and (2.9), we see that we have two methods of obtaining the 
transfer function of our filter. Since the transfer function completely characterizes the 
behavior of the filter and, hence, the ocean medium that the filter is modeling, we need 
some method of determining this quantity. Equation (2.9) is not much use for our 
problem since it assumes that we already know the impulse response of the filter which we 
do not. However, Eq. (2.8) is much more useful. It says that if we apply a plane wave 
input to the system, then there is a simple relationship between the input and output that 
determines the value of the transfer function. However, this relationship is more versatile 
than may first appear. It tells us that if we can express an output due to any type of input 
as the product of the time-harmonic plane wave given by Eq. (2.6) and some other factor, 


then this other factor is the transfer function. 


Expressing the impulse response term in Eq. (2.5) in terms of the two dimensional 
inverse Fourier transform of its transfer function and combining with the input term, it can 
be shown that the two dimensional Fourier transform of the output of a linear time- . 


invariant, space-variant filter, given by Eq. (2.5), has the form 


Y(f,B) = f f X(f,v)H(E,v;r) 12" o ay ar (2.13) 


which is the frequency and angular spectrum of the output. Note that the output spatial 
frequency term is different, in general, from that at the input, and that Y#XH. The 
exception to this rule is when the ocean medium is homogeneous. In this case the filter 
model is space-invariant and Eq. (2.13) collapses to the form Y= XH , which is the same 
as the familiar result encountered for one-dimensional LTI systems. 

We can determine the impulse response once we know the transfer fountain by using 


the inverse transform given by 


A(t gsr) = f f H(f.vsr) 7 @ Po at av , (2.14) 


which can be expressed in transform notation as 


h(t,rg3r) = Fp Fy ( H(fv:r) } (2.15) 


where Fy represents the inverse time-domain Fourier transform and lak represents the 
: é : - -1. 
inverse spatial Fourier transform. The form of Fr and Fy is the same as for Egs. (2.11) 


and (2.12) with the sign of the exponent changed and the variable of integration changed to 
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f and v, respectively. This is the same as for the familiar one-dimensioned time-frequency 
Fourier transform. 

Note that the development in this section has so far been completely general and can be 
applied to any linear, time-invariant, space-variant system. We have not made any 
assumptions regarding the space variance of the system, i.e., we have not specified a 


particular form for the sound speed-profile as a function of three-dimensional space. 


C. THE PROBLEM STATEMENT 

Let us now consider the underwater acoustic problem from the viewpoint of linear 
systems theory. Our aim is to be able to compute the output electrical signal from a receive 
transducer due to an applied electrical signal at a transmit transducer. Simply stated, this 


process consists of three steps : 


* conversion of an electrical signal into a transmitted acoustic signal, 
* propagation of an acoustic signal through the ocean medium, and 


* conversion of a received acoustic signal into an electrical signal 


Our intention is to model each of these steps as a time-invariant, space-variant / invariant 
filter as appropriate. The system model conforming to this plan is shown in Fig. 2.2 in 
block diagram form. 

The remainder of this chapter will deal with the development of expressions describing 
the behavior of the transmit and receive arrays, and the relationship or coupling equation 
that relates the output, y(tr), to the input, x(t,r). The development of the transfer function 


that characterizes the ocean medium will be discussed in the next chapter. 






ym (tr) 






hy (7,rg;r) 
Hy (f,v;r) 


transmit ocean receive 
array medium array 


Fig. 2.2. System Model 


D. COUPLING THE TRANSMITTED SIGNAL TO THE MEDIUM 

Consider now the specific problem of transmitting an electrical signal through the ocean 
medium. The electrical signal will be used to drive some transducer (or array of 
transducers) which will convert the electrical energy into acoustic energy. We therefore 
need to find a relationship that describes this transformation of an electrical signal x(t,r) 
into an acoustic signal xjy(t,r). The acoustic signal represents the rate at which fluid 
volume is added at time t and position r per unit volume of fluid. 

Suppose we model an infinitesimal region of some arbitrarily shaped transducer as a 
linear filter with a corresponding impulse response. The output from such a transducer at 
a given position r would be given by the time-domain convolution of the corresponding 
impulse response with the input signal. Translated into the frequency domain via a Fourier 


transform with respect to time, this relationship can be expressed as 


XM (f,r) = X(F)AT (f,r) ; (2.16) 


where Ay (f,r) is the complex frequency response (or aperture function) of the transducer 


at spatial location r. Taking the spatial Fourier transform of both sides of Eq. (2.16), it 


can be shown that: 


Xm (fv) -f X(f,a)D 7 (f,v-a) da (2.17) 


=Se 


where Dy(f,a.) is the spatial Fourier transform of the complex frequency response of the 
transducer and is known as the far-field directivity function or beam pattern. Note that the 
result given by Eq. (2.17) is valid for a completely arbitrary volume aperture. 

If we make the simplifying assumption that an identical electrical signal is applied at all 


spatial locations, r, of the transducer then 
M(ET) = x(0): (2.18) 


Taking both the time and spatial domain Fourier transforms of a signal of the form of 


Eq.(2.18) and substituting into Eq. (2.17) yields 
Xv (f,v) = X(D (Fv). (2.19) 


If the transmit aperture is a planar array, then for most practical situations we can consider 
that an Senne input electrical signal is applied to all the elements of the array before any 
electronics used for beamsteering. As a result, the restriction imposed by Eq. (2.18) does 
not impose any limitations. What we gain, however, when comparing Eqs. (2.17) and 


(2.19), is a significant simplification in our analysis. [Ref. 4:p. 4] 


E. THE TRANSMIT ARRAY 
Assume that the transmit aperture is a planar array of identical transducer elements. 


The product theorem for planar arrays, given by 
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D(f,v) = S(V)E(E.V), (2.20) 


states that the far-field directivity function of such an array is equal to the far-field 
directivity function of one of the transducer elements, E(f,v) , multiplied by the far-field 
directivity function of an equivalent array of complex weighted point sources, S(v). We 
will assume for the purposes of this thesis that our transmit array is located in the XY plane 
and consists of MT x NT (odd) complex-weighted point sources. We will further assume 
that the array elements are equally spaced in the X and Y directions; however, it should be 
noted that this is not a constraint imposed by the use of the product theorem. 

Since the spatial Fourier transforms deal with position relative to some source, it 
makes sense to locate the transmit array in a convenient position to simplify the analysis. 
Accordingly, we position the array such that its center forms the origin of the coordinate 
frame of reference in the X and Z directions while the ocean surface forms the center in the 
Y direction, i.e., the array is centered at (x, =0,y9 =yrT.Z) =0). The geometrical 
interpretation of this configuration is shown in Fig. 2.3, which is also the physical 
representation of the block diagram of Fig. 2.2. 

Given these conditions, Ziomek [Ref. 3:p. 126] shows that the far-field directivity 
function of the transmit array is given by 


MT" NT 
j2nfy md j2nfynd j2nf 
Dr(fvy= YY calf XM XT MYMYT QM YT (2.21) 


where Cm_(f) is the frequency dependent complex weight at the transmit element (m,n), 


and dy and dyy are the interelement spacings in the X and Y directions, respectively. 
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Fig. 2.3 Transmit and Receive Array Geometry 


The last complex exponential is a phase factor that accounts for the array not being centered 


at the origin. Also, 


MT =(MT-1)/2 ( 2.22a) 


and 


NT =(NT-1)/2. (2.22b) 


It should be noted from Eq. (2.21) that a single omnidirectional point source and a linear 


transmit array are special cases of the planar transmit array. This is also what we would 


expect from examining this problem from a purely physical point of view. If we assume 


that the transmitted electrical signal is applied to the array elements before the complex 
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weighting, then we can say that an identical electrical signal is applied to all the elements in 
the array. Accordingly, Eq. (2.21) can be combined with Eq. (2.19) to describe the 
acoustic input to the ocean medium. 
Finally, what is meant by the far-field directivity function? It can be shown that a field 
point at a range r is considered to be in the far field if 
mR? 


where R is the maximum radial extent of the aperture. For our planar transmit array, this 


distance is given by 
5 1/2 
R=[(MT xp) +(NTdyp) 1] (2.24) 


The physical significance of being in the far field is that the range is large compared to a 
wavelength and the separation between the transducer elements. Consequently, any 
second or higher order terms in the expression for the directivity function can be ignored 


and the mathematics of the problem is simplified. [Ref. 3:p. 38] [Ref. 2:p. 170] 


F. THE PHASED ARRAY 
Now let us take a closer look at the frequency-dependent complex weighting function 


Cmn (f) of the far-field directivity function given by Eq. (2.21). Since this weighting 


function is complex, it can be expressed as 


Cmn (f) = Amn (f) e!@mn®) (2.25) 
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which consists of an amplitude weighting component am, (f) and a phase weighting term 


©mn(f). Now consider that this phase weighting is a linear phase variation of the form 
Omn (f) = -20 (fy mdyz + fy ndyy) (2.26) 


where m, n, dxz, and dyr are defined as for Eq. (2.21). If we define the spatial 


frequency terms in Eq. (2.26) as 


fa (2.27a) 


and 
ey cites (2.27b) 


then, substitute Eqs. (2.25) and (2.26) into Eq. (2.21), it can be shown that the resulting 
far-field directivity function has been steered in the direction u = u' and v = v'.[Ref. 3:pp. 
132-134]. 

The term phased array means nothing more than an array which employs beamsteering 
using phase weighting. Inspection of our assumed form of the phase weighting term 
shows that it is both a separable function and an odd function of the indices m and n. 
These are necessary conditions in order to carry out the beamsteering. It should be noted 
from the preceding development that no assumption was made as to whether the amplitude 
term was either separable or an odd function of the indices m and n. For the purposes of 
this thesis, we will only be considering rectangular amplitude weighting, which is a 


separable function and leads to a simplification in the analysis. 


17 


Since the phase weighting is already a separable function, the assumption that the 
amplitude weighting is also separable results in the entire complex weighting term given by 
Eq. (2.25) being separable. Substitution of this form of weighting into Eq. (2.21) results 


in the expression for the far-field directivity function being separable which gives 
Dr (f,v) = Dr (ffx) Dr (ify) (2.28) 
where 


MT' : : 
Dr(ffx)= >, Ae?@xmixr . ity maxT (2.29) 
m= -MT' 


A is the magnitude of the amplitude weighting, and a similar expression exists for 


Dr(f,fy). To give a more intuitive feel for what this far-field directivity function looks 
like, we can express Eq. (2.29) in the alternative form 


sin[n(fx - fy MT xq] 


Dr(ffy) =A (2.30) 


sinfn(fy - fy dq] 


which looks similar to a sinc function when viewed in direction cosine space. The 


expression for Dy (f,fy ) is of the same form as Eq. (2.30) except that it is multiplied by the 
last complex exponential phase factor in Eq. (2.21) that accounts for the array not being 


centered at the origin. 
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G. RECEIVED ACOUSTIC FIELD 
We next need to calculate the acoustic signal at some spatial location r = (x,y,z) after 
transmission through the ocean medium. Applying the space-variant, time-invariant 


impulse response of the ocean medium to Eq. (2.5) yields 


yen = f if Xm(t-t.r-Po hy (t.rosrdtdro . @:31) 


oO —oo 


If we express the input signal in the previous expression in terms of its frequency and 
angular spectrum and collect terms which combine with the impulse response to produce 
the transfer function of the ocean medium, we obtain 


-j2uver 


Yu (f,r) = J XM (f,v)Hy (f,v;r)e dv (2.32) 


oOo 


which is the frequency spectrum of the output acoustic field and is one of our desired 
results. 

Next we need to consider some receive aperture to convert the acoustic energy into 
electrical energy. For the purposes of this thesis, we will assume that the receive aperture, 
like the transmitter, is a planar array that is characterized by a complex aperture function 
Ag (f,r). Analogous to the development of Eq.(2.16), it can be shown that the 


corresponding output electrical signal is given by 


y(t,r) = le Yu (BnAgiine df. (2.33) 
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We will assume that our receive array is located in a plane parallel to the XY plane and 
consists of MR x NR (odd ) point receivers centered at (XR, YR>YR). We will further 
assume that the array elements are equally spaced in the X and Y directions. Due to the 


assumption of point receivers, the aperture function can be expressed as 


MR’ NR’ 
Agir)= > & wmnd[x-(xp tmd xp ISLy-(yp tndyp 8[z-zR] (2.34) 
m=-MR’ n=-NR’ 


where W mn (f) is the frequency dependent complex weight at element (m,n), dxp anddyr 


are the interelement spacings in the X and Y directions, respectively, 


MR'=(MR-1)/2, (2.35a) 


and 


NR'=(NR-1)/2. (2.35b) 


For this thesis, we are not concerned with processing the received electrical signals. 
What we want to accomplish is the ability to compute the complex acoustic field at the 
receive array as a function of frequency and spatial coordinates, Eq. (2.32), and the time- 
domain electrical signal at each element in the receive array. The signal processing that 
will be done on the computer simulated received electrical signals will be via auxiliary 
programs, e.g., FFT and adaptive beamformers. Accordingly, we are primarily concerned 


with formatting our outputs in a form acceptable to those programs. 
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As a consequence of these requirements, and in an attempt to minimize computation 
time, we will assume the simplest form of the aperture function where the complex 
weighting function W mp (f) is unity. Applying this simplifying assumption to Eq. (2.34) 


and substituting the resulting form of the aperture function into Eg. (2.33) yields 


y(t,m,n,zp ) = Yy (f,m,n,zp er df (2.36) 
R M R 


which is the time domain electrical signal at element (m,n) and is our other desired result. 


H. OVERALL SYSTEM COMPLEX FREQUENCY RESPONSE 
If we combine Eqs. (2.19) and (2.32) we obtain 


Ym (fr) = X(MH(f,r) (2.37) 
where 


H(fr) = f Dr(f,v) Hy (Evin av (2.38) 


is known as the overall system complex frequency response [Ref. 4:p. 5] [Ref. 6:p. 1672]. 
Note that in the development carried out in this chapter we have not made any assumptions 
about the ocean medium transfer function. The only assumption used in the derivation of 
Eq. (2.35) was that an identical electrical signal was applied at all spatial locations of the 


transmit array before the application of the complex weights. 


2M 


Taking the inverse time-domain Fourier transform of Eq. (2.37), and using the 


expansion of Eq. (2.38) it can be shown that 


ym (tr) = i" if X(ND (fv) Hy (Evie 7" avd (2.39) 


which is the coupling equation that relates the output acoustic signal from the medium, 


YM (tr) , to the frequency spectrum of the applied elctrical signal X(f). 
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I. THE AXISYMMETRIC CASE 


A. OVERALL SYSTEM COMPLEX FREQUENCY RESPONSE REVISITED 
Let us examine the expression for the overall system complex frequency response, 
given by Eq. (2.38), that was developed in the previous chapter. It would appear from an 
initial inspection that Eq. (2.38) is a three-dimensional integration with respect to 
dv =dfydfydfz . However, it should be noted that since the integrand is a function of 
the frequency, f, as well as the spatial frequencies, v = (fx,fy,fz), the integral can be 
reduced from a three-dimensional to a two-dimensional integral. This is due to the manner 
in which the spatial frequencies given by Eq. (2.7) are related to the frequency f via 
direction cosines. The property of direction cosines that the sum of their squares always 
equals unity means that, since we know f, we can express one of the direction cosines in 
terms of the other two. Thus, we do not need to integrate with respect to all three spatial 
frequencies in order to span direction cosine space. Using the spatial frequency definitions 


of Eq. (2.7), it can be shown that 


H(f,x,y,z) = f i Dr(f.fx .fy .fz)HM (ffx sly iz s%.Y.Z) 


oOo = =0O 


x SIE RE QP fyy + £72) dfx dfz (3.1) 


where 


1/2 
fy =£[(f/c9)-(fx +f2)] . (FR +f2) S (flog) (3.2a) 


and 
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1/2 
fy = Fi1(& +f) (elegy) Cheeta eee (3.2b) 


and Co is the speed of sound at a particular transmit array point source element. The plus 
(minus) sign in Eq. (3.2a) is chosen whenever y 2 yo (y < yo) and corresponds to the 
field point being deeper (shallower) than the source point, i1.e., downward (upward) 
traveling waves. The minus (plus) sign in Eq. (3.2b) corresponds to the plus (minus) 
sign of Eq. (3.2a) in order to generate evanescent waves. 

Note that we have not made any assumptions about the ocean medium transfer function 
in our development up until this point. However, for the purposes of this thesis, we are 
going to assume that the acoustic energy of the transmitted signal is contained within a 
sound channel, e.g., the SOFAR or deep sound channel. In this situation the sound 
propagates between upper and lower turning points with no surface and bottom boundary 
interaction. We therefore expect the acoustic field trapped within this channel to exhibit 
some form of cylindrical spreading. If this assumption is correct, then cylindrical 


coordinates should provide an easy reference frame against which to describe this behavior. 


B. THE CYLINDRICAL COORDINATE SYSTEM 
The cylindrical coordinate system used in this thesis is shown in Fig. 3.1, where the 


transformation equations that relate the rectangular and cylindrical coordinates are given by 
x =rcosd, y=y,andz=rsino (3.3) 
where r =V x? + 2 is the polar radial distance. 


In addition, the spatial frequencies must be transformed using the following set of 


equations ( see Fig. 3.2 ): 
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xX 


—_ if = (1,,y) 
(cross range ) 


ry. 
(depth ) 


Fig. 3.1 The Cylindrical Coordinate System. 





a6 


Fig 3.2 3-D Wave Propagation. 


2 


fx =f,cos, fy =fy, and fz =f,sinG (3.4) 
where f, = V f a +f ; is the spatial frequency in the radial direction. As aresult, 
dfx df7z > f,df,d. (375) 
The physical significance of the two different angles in Eqs. (3.3) and (3.5) is that relates 
to the position of some field point whereas € relates to the direction of propagation, and 
they are not necessarily the same. 


Applying these definitions to Eq. (3.1), it can be shown that the overall system 


complex frequency response can be expressed in cylindrical coordinates as 


2 oo 
H(f,r) os ii ii Drtt; siete )Hy (i faGaiy, sr) 
0 0 


xe @ 2nllar cose sing + fyy + fer sing sind) + ge ge ae 
where 
fy ee) ee  aGuie.)). (3.7) 
and 
fy = Fil f “(f/eo)1 f >(f/eg)*, (3.7b) 


and Cg is the speed of sound as previously defined. 
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Now consider the specific transmit array geometry discussed in Chapter 2. If we 
substitute the far-field directivity function for this array, as given by Eq. (2.21), into the 
previous equation, and perform the necessary transformation from rectangular to 


cylindrical coordinates, we find that 


2n oe MT’ NT 
H(f,r) = f ih » SY cma) Hy (,f,.6,fy sr) 
0 0 m=-MT n=-NT 


x ¢ HM r(t cost sing + r sing sing + mdx poosh) | -j2nfyly - yp - nyt) ¢ df dC . G8) 


C. THE AXISYMMETRIC OCEAN MEDIUM TRANSFER FUNCTION 

In this thesis, we shall assume that the speed of sound in the ocean medium is a 
function of depth only. It will be shown in the following chapter that under this condition 
the ocean medium transfer function is axisymmetric, i.e., independent of the azimuthal 


angle € and, as a result, it can be expressed in the form 
Hy (ff, .¢,fysr.0.y) = Hy (ff, .¥93y) (3.9) 


where yg =y7 - ndyz is the depth of a point source in the nth horizontal row of the 
transmit array. Note the presence of the terms y and yo in Eq. (3.9). These terms 
indicate that the ocean medium transfer function is dependent on the speed of sound at the 
position of the transmitter as well as the receiver. Since the assumed form of our 
transmitter is a planar array, we need to consider the depth and, hence, speed of sound at 
each element in the array when calculating the acoustic field at some spatial location r. This 
makes physical sense when considering Eq. (3.2), which shows that the speed of sound at 


the transmit element affects the direction of propagation and the threshold between 


e]/ 


propagating and evanescent waves. If the form of the ocean medium transfer function 


given by Eq. (3.9) is substituted into Eq. (3.8), then the integrals can be regrouped to give 


i -j2nfyly - 
H(f.r) = i! 3 oe Conf) Hy Eofpsyosy) Cem YY Yo ¢, 
0 m=-MT n=- 


fia eee cos sind + rsin€ sind + md x TOD) ar gg. (3.10) 
0 


Next, consider the integration with respect to in the preceding equation. Due to 
our placement of the transmit array in the coordinate system, the angle @ is the horizontal 
polar angle (relative to the positive X axis) between the center of the transmit array and the 
field point of interest. Now define a new angle © » that is the horizontal polar angle 
between a specific array element (m,n) and the field point, as shown in Fig. 3.3. Note that 
the field point in Fig. 3.3 is in three-dimensional space, but that it is only the (x,z) 
component of that position that determines 6. Also note that since the angle is in the XZ 
plane, it is independent of the vertical position y and, hence, is independent of the array 
index n. Expressing the horizontal polar angle 9 in this integral in terms of the new angle 


%m and collecting terms, results in 


an, : Pe an, 
f inte cosC sing + r sin sing + mdy ;cosl) at - LY oe J2mf tr mcos( - Om) ge (3.11) 
0 0 
where 
2 @ Me 
Im =[(K-mdxy) +z ] (3.12) 
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Fig 3.3. Graphical Interpretation of > ,, ( Plan View ). 


is the polar radial distance between element (m,n) in the transmit array and the field point, 


and md xr is the horizontal distance between element (m,n) and the center of the array. 
Brekhovskikh [Ref. 8:p. 76] and Officer [Ref. 9:p. 126] show that the right-hand side 
of Eq. (3.11) is in the form of the integral representation of the zero-order Bessel function 


of the first kind where 


75 Se 
if @ 12m Tmoos( : om ae = 2nJ g(2nf tm) - (3.13) 
0 


Substituting Eqs. (3.11) and (3.13) into Eq. (3.10) yields 


co MT" ae 
Hic) = 20 i > D  Cmnlf) He (ff ¥osy) Jo(2nf rtm) 
0 


x eM -Ye ge (3.14) 


pe) 


which is the fundamental result of this chapter. Comparing the previous equation with Eq. 
(3.1), it can be seen that our expression for the overall system complex frequency response 
has been reduced from a double to a single integral. This was based on the assumption 
that the speed of sound was a function of depth only and, as a result, the corresponding 
ocean medium transfer function was axisymmetric. Since numerical integration is a time 
consuming process, it is surmised that the reduction in the number of integrations required 
in the computation of our computer-generated acoustic output should result in the saving of 
a significant amount of computer time. 

It should be emphasized that Eq. (3.14) is valid for any axisymmetric ocean medium 
transfer function. Since this condition is satisfied for the case when the sound speed is a 
function of depth only, Eq. (3.14) represents the general form of the overall system 
complex frequency response that we will use in the development of the next chapter. 
Remember that we are trying to develop a general, modular, pulse-propagation model 
based on sound-speed profiles that are a function of depth. Thus Eq. (3.14), when 
combined with the coupling equations of Chapter 2, is the general solution for this 
assumed type of sound-speed profile. In the next chapter we will develop a specific 


solution based on the WKB approximation. 


D. AXIAL SYMMETRY OF THE TRANSMIT ARRAY 

The only assumption made regarding axial symmetry in this development was that the 
ocean medium was axisymmetric. This makes physical sense when considering the ocean 
medium itself, in the absence of any source of acoustic signal. Since the speed of sound 
was assumed to be a function of depth only, there is nothing in the medium that should 
cause any azimuthal } dependence on the propagation of an acoustic signal. Accordingly, 


the ocean medium transfer function should be independent of the azimuthal angle 9. 
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The assumption of axial symmetry may at first glance appear to be invalid when we 
include the transmit array in the problem, as the beam pattern or far-field directivity 
function of the planar array is not generally symmetrical about the Y axis. However, 
consider the problem from the following perspective. From linear systems theory, the far- 
field directivity function of the array is the linear superposition of the far-field directivity 
functions of the individual array elements. The far-field directivity function of a point 
source is unity, i.e., itis omnidirectional. The far-field directivity function of a complex 
weighted point source is also omnidirectional, since the complex weight factor is merely a 
scaling term. Equations (2.20), (2.21) and (2.29) demonstrate that the far-field directivity 
function of the transmit array is a result of the way these omnidirectional beam patterns are 
combined rather than by altering the beam pattern of the individual point sources. The 
effect of the complex weighting is to control the way these individual far-field directivity 
functions are combined. 

Examination of Eq. (2.21) also shows that there is an additional phase term that is 
associated with each array element that accounts for that element's displacement from the 
origin of the coordinate axis system. Thus, we can assume that the far-field directivity 
function of the transmit array is the linear combination of complex weighted axisymmetric 
point sources. Accordingly, the overall system complex frequency response solution 
given by Eq. (3.14) can be thought of as applying to each individual array element, with 
the solution for the array being a linear combination of these results. This is nothing more 
than rewriting Eq. (3.14) with the summations outside of the integral. Since the 
integration and summations are with respect to different quantities, regardless of whether 
the integration is inside the summations or outside, the value of the overall system complex 


frequency response will be the same. 
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IV. THE OCEAN MEDIUM TRANSFER FUNCTION 


A. THE HELMHOLTZ WAVE EQUATION 

Examining the various coupling equations developed in Chapter 2, we see that they 
already incorporate an input acoustic signal x,y (t,r), with corresponding frequency and 
angular spectrum X yy (f,v). This means that the transfer function that characterizes the 
ocean medium, Hy (f,v;r), is independent of the input as we would expect from linear 
systems theory. Accordingly, when trying to find a transfer function to describe the 
ocean medium we only need to solve the wave equation given by 


22 
Vo(t.r) - ae o(t.r) = 0 (4.1) 
c (r) ot 


instead of the linear acoustic wave equation given by Eq. (2.1). 
The first assumption that we will make in finding a solution to Eq. (4.1) is that the 
velocity potential, @(t,r), has a time-harmonic dependence, that is, 


j2nft 


(tr) = P(re (4.2) 


If this form of the velocity potential is then substituted into Eq. (4.1) and the time-harmonic 
dependence is factored out, the result is the time-independent Helmholtz wave equation 


given by 


V‘o(r) +k2(ro(r) = 0 (4.3) 
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where k(r) = 2nf/c(r) is the wave number and c(r) is the speed of sound in the fluid 
medium at spatial location r. 

As discussed in the previous chapter, we will be assuming for the purposes of this 
thesis that the speed of sound is a function of only one spatial variable, the depth y, i.e., 
c(r) = c(y). While this assumption is not in general true, it is nevertheless adequate 
enough to describe the ocean medium for most circumstances where we are not considering 
propagation ranges longer than say, 20 km, or propagation across significant 
oceanographic perturbations such as fronts or eddies. We can also consider that the 
sound-speed profile in the region of the deep sound channel is less variable than for 
shallow sound channels or surface ducts. Thus, the assumption that the sound-speed 
profile is not dependent on geographical position is adequate for our model, within the 
limitations discussed. 

This assumption allows us to make the substitution k(r) = k(y) in Eq. (4.3) which 
greatly simplifies the analysis. Since k(r) is now a function of only one variable, we can 
use the method of separation of variables to solve for the acoustic velocity potential. 
Applying this solution technique in cylindrical coordinates, we assume a solution for @(r) 


of the form 
p(r) = RDO) Y(y) (4.4) 


where r=(r,@,y) is the spatial location in cylindrical coordinates (see Fig. 3.1) and the 
Laplacian appearing in Eqs. (4.1) and (4.3) is defined as being 
oo at a ee 
——— 
G] 2 
r or 


; (4.5) 
rao dy 
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B. AZIMUTHAL ANGLE AND RANGE TERM SOLUTION 

Since the speed of sound is a function of depth only, there is nothing in the medium 
that should cause any azimuthal, i.e., 6, dependence on the propagation of an acoustic . 
signal. Accordingly, the solution for the acoustic velocity potential should be independent 
of the azimuthal angle @. This is referred to as the ocean medium having axial symmetry. 


Applying this assumption, Eq. (4.4) reduces to 
p(r) = R@Y(y) (4.6) 


and the partial differential with respect to azimuthal angle in the Laplacian of Eq. (4.5) 
drops out. 

If we substitute Eq. (4.4) into Eq. (4.3) and assume that k(r) = k(y), the method of 
separation of variables yields 


a 
ae aoe +k,R(r) =0 (4.7) 
or 5 


where the wave number in the radial direction, k,, has been used as a separation constant. 
Making the substitution R(r) = g( k,r ) into the previous equation, it can be shown that Eq. 
(4.7) reduces to Bessel's differential equation of order zero. The solution to this equation 


can be expressed in either of the forms 
R(r) = AJg (k,r) + B Yo (k,r) (4.8a) 


or 
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R(t) = AH} k,1) + BHO Xk, 1) (4.8b) 


where Jg and Yg are the zero order Bessel functions of the first and second kind, 
respectively, with Yo also being known as the zero order Neumann function. The 
functions HY? and He are zero order Hankel functions of the first and second kind, 
respectively, or Bessel functions of the third kind. Since Eqs. (4.8a) and (4.8b) are 
alternative expressions describing the same quantity, the Hankel functions should be able to 
be expressed in terms of the Bessel functions of the first and second kind and visa-versa. 
This is, in fact, the case and we will make use of one of these relationships at a later stage 
of the analysis. 

Since there are two alternative representations of the range term R(r), we would expect 
that one representation is more useful than the other for the situation that we wish to model. 


If we consider the approximations for large argument for the two Hankel functions 


1 / 2) i(k, -n/4 
Hi 1k, i) =e re el ay (4.9a) 
I 


and 


2 i(k, m4 
HS k,n) > 4f ees (4.9b) 
8 


then because of our assumed form of the time-harmonic dependence of the velocity 
potential given by Eq. (4.2), Eq. (4.9a) represents an incoming wave while Eq. (4.9b) 
represents an outgoing wave. Since we are assuming that our transmitted acoustic signal 


is traveling in a channel with no boundary interaction, we will have only outgoing waves 


Se) 


from the source and no incoming waves due to reflection. It then appears that Eq. (4.8b) 
with the constant A set to zero, is an appropriate solution for our range term. We will 


therefore assume that our range term solution is of the form 
~ (2) 4 
R(r) =B,Hg’ (k;r) (4.10) 


where B, is a constant. By comparison, if we look at the approximations for large 
argument of the Bessel functions of the first and second kind, we find that they are 
comprised of both incoming and outgoing wave components and so are not well suited to 
our problem. 
Upon substitution of Eq. (4.10) into Eq. (4.7) however, we find that our assumed 
form of the solution for the range term is a solution everywhere except at the origin, i.e., 
(2 


r= 0, since H¢ ) “blows up" at this point and equals -joo. Using a volume integral 


approach, Gauss's theorem, and the asymptotic expression for the Hankel function given 





by 
Jim | HO «,1) == In (k,r), (4.11) 
it can be shown that Eq. (4.10) is a solution to 
ioe kaR Ge ioe Sir) (4.12) 


for all r, where Eq. (4.12) is the same as Eq. (4.7) with the addition of a forcing function 


or source term. This makes physical sense when we remember that Eq. (4.10) describes 
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outgoing waves in the radial direction. Since these outgoing waves have to come from 
somewhere and there are no incoming waves, there must be a source of some description, 
that generates the waves. Thus Eq. (4.10) is a solution to Eq. (4.12) that satisfies the 
boundary condition of a point source at the origin. The choice of a point source as the 
forcing function is appropnate since we have assumed that our transmit aperture is a planar 


array of point sources. 


C. THE WKB APPROXIMATION 

For the depth dependent term Y(y) of our solution for the velocity potential, we will 
obtain an approximate solution based on the geometrical optics approximation, otherwise 
known as the WKB (Wentzel, Kramers, and Brillouin) approximation. If we substitute 
Eq. (4.4) into Eq. (4.3) and assume that k(r) = k(y), the method of separation of variables 
yields 


2 
[ 4+ | vor=0 (4.13) 
oy 


where the wave number in the vertical (or depth) direction, ky (y), is calculated using the 


relationship 
Di 
Ke (yy = ky) -k? (4.14) 


If the vertical wave number was a constant, i.e., ky (y) = ky, meaning that it was no 


longer a function of depth, then Eq. (4.13) would have the exact solution 


Y(y) = Ae T*Y" 4 Belly”, (4.15) 


iy 


If we now assume that the local medium is homogeneous, that is the variations of the 
properties of the medium per wavelength are small, then a solution similar in form to Eq. 
(4.15) would be a good approximation. The WKB method assumes that the solution to 
Eq. (4.13) is of the form 


Y(y) = Ay (ye (4.16) 


where A y (y) and @y (y) are real amplitude and phase functions respectively. Substituting 


this assumed form of our solution into Eq. (4.13), and making the simplifying assumption 


that 


JAY) /Ay@)|«lO¥Q) 1" (4.17) 


yields 


y y 
k j k 
I roa, if yO) at Jacacnt 


il a 
*O)= Tky Gt [ Ave Yo 


which is the approximate WKB solution [Ref.3:pp. 208-220] [Ref. 8:pp. 129-134]. It 
should be noted that the only assumption used that leads to this not being an exact solution 
is that given by Eq. (4.17). This method also assumes that there is no reflection in a 
horizontally stratified medium. 

The LHS of Eq. (4.17) is the relative differential of the time rate of change of 
amplitude and is a measure of the rate of change in the properties of the medium. If the 


medium was homogeneous, then this term would be zero since the characteristics of the 
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medium would not be varying. Therefore this term is a measure of the inhomogeneity of 
the medium. The RHS of Eq. (4.17) is inversely proportional to wavelength since for a 
given frequency the rate of change of phase increases with decreasing wavelength. 
Equation (4.17) is therefore a reasonable simplifying assumption since it is nothing more 
than a mathematical formulation of the assumption made earlier that the variations of the 
properties of the medium are small compared to a wavelength. 

The integrals in Eq. (4.18) give the phase change as the waves propagate between 
depths y and yo, while the factor in front is to preserve the conservation of energy. Due to 
the assumed form of the time-harmonic dependence of the velocity potential given in Eq. 
(4.2), the first and second terms in Eq. (4.18) correspond to waves traveling in the positive 
(down) and negative (up) Y directions, respectively. It should be noted that the WKB 
solution is then the superposition of two waves traveling without interaction in opposite 
directions. Also note that the approximate solution is not valid as ky (y) approaches zero 
since the solution approaches infinity. The point at which this occurs is referred to as a 


tuming point. 


D. THE VELOCITY POTENTIAL SOLUTION 
Combining Eqs. (4.6), (4.10) and (4.18), we obtain a solution for the velocity potential 


of the form 
AHMan f key at 
ee (4.19) 
ky (y)| 
where 
oe aa 2 2 
ky (y) = +2n[ (f/cly))- fy] wCi e SGC), (4.20a) 


39 


and 


2 2 ii 2. 2 
ky (y) = ¥j2a[f,-(f/cty)) ] > C§PPHCl/cOyr, (4.20b) - 


which is our separable function solution to the wave equation. The constant A in Eq. 
(4.19) combines the constants Ay and B,; of the preceding equations. Note that we have 
replaced the approximation symbol with the equality sign in Eq. (4.19). This is only a 
matter of notation, and we still recognize that the solution is only an approximation. 

Since the linear wave equation is in terms of partial differentials with respect to y and r, 
if we multiply Eq. (4.19) by a function of k,, then the result is also a solution to the wave 
equation since k, is independent of y and r [Ref. 9:p. 125]. Similarly, any summation of 
such solutions will also be a solution to the wave equation. Combining these concepts into 


a generalized integral form results in 


_ (2) : ii - 
oa f Ak) HO) JJ kKyO oo 


a EC 
aky (y)| ie ; 


where A(k,) is the term that incorporates the constant A of Eq. (4.19) and the arbitrary 
function of k, as discussed. Equation (4.21) is therefore our solution of the linear wave 
equation based on the WKB solution for an inhomogeneous medium. If Eq. (4.21) is 
correct, then for the case of a homogeneous medium it should collapse into the free-space 


Green's function for an isospeed medium. As well as validating the form of the equation, 


it should also provide us with an expression for the term A(k, ). 
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E. THE FREE-SPACE GREEN'S FUNCTION 

The free-space problem for an isospeed medium is one in which there are no 
boundaries to consider and the medium is considered to be homogeneous. Since one of 
the assumptions made in the development of Eq. (4.21) was that there was no surface or 
bottom boundary interaction, all that we have to do to find a solution for the free-space, 
isospeed problem is to apply the conditions for a homogeneous medium to Eq. (4.21). 
For the homogeneous medium, ky (y) = ky(y,) is a constant. For this case, there is a 
closed form solution of the integral appearing in the exponent of Eq. (4.21), which results 


in a solution of the form 


g(r) = f Aer) HOCkr2) ky Gg DEN ee (4.22) 


— AkyGo)) 


It should be noted that for the case of a homogeneous medium, the simplifying assumption 
given by Eq. (4.17) is unnecessary, and Eq. (4.22) is an exact solution. 

We will now specifically solve the time-independent Helmholtz equation of the 
form of Eq. (4.3) for the free-space, isospeed problem and compare the result with Eq. 
(4.22). In the free space problem a wave will radiate in an outwards direction, 
experiencing spherical spreading, and be free of boundary interactions. The homogeneous 
medium condition means that k(r) = k is a constant. Since the wave will spread 
spherically, it makes sense to use the spherical coordinate system to describe the solution. 


It can be shown that 





Or) = (4.23) 


4] 


is a solution to Eq. (4.3), under the conditions described, everywhere except at r = 0. 
Note that in Eq. (4.23) r is a scalar quantity rather than a vector. Equation (4.23), 


however, does satisfy 
2 2 
V o(r) +k g(r) = -4xA&(r) (4.24) 


for allr. This is consistent with there having to be a source at r = 0 to generate the wave. 
The choice of a point source as the source is a convenient one since in Chapter 2 we 
assumed that our transmit aperture was an array of point sources. If we think of this 
forcing function or source term as having unit amplitude, we set the constant in Eq. (4.23) 
equal to (-1/4n). Applying this amplitude term and converting to cylindrical coordinates 


we can show that 


ps 
A a ana (4.25) 
is a solution to 
V “o(r) +k “o(r) = oe Se B(y- -Yo) (4.26) 
for all r where 
2 Me 
R={t +7 Vo)) ee: (4.27) 
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and r is the polar radial distance. Equations (4.23) and (4.25) are the free-space Green's 
functions for an isospeed medium in spherical and cylindrical coordinates, respectively. 
They represent the solution to the time-independent Helmholtz equation for a homogeneous 
medium with no boundary interactions, and matching the boundary condition for a point 
source. 

Officer [Ref. 9:p. 126] shows that the free-space Green's function in cylindrical 
coordinates can be expressed as an integral of the zero-order Bessel function of the first 
kind, having the form 


-jkR oy 
e & ie jkyVo) (Y-Yo) 
R J Jky (Yo) J (kyr) € dk; . (4.28) 


If we now combine Eqs. (4.28) and (4.25) we obtain 


ik; iky Vo) -Yo) 
Pry S Jo(k,ne 2 o’ dk 4.29 
which is the solution found by solving the linear wave equation specifically for the free- 


space, isospeed problem. 


F. DETERMINATION OF THE OCEAN MEDIUM TRANSFER FUNCTION 
Returning to our solution for the free-space problem based on the WKB approximation 
given by Eq. (4.22), we see that it is in a similar form to Eq. (4.29) but is in terms of a 
Hankel function rather than a Bessel function. Based on observations made earlier in this 
chapter, however, we would expect to be able to express Eq. (4.22) in terms of a Bessel 


function. This is in fact the case, and using the relationships Ho (-k,1) = -H (k,1) and 
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Jo (kyr) = 5 Lo kt) + HO k,n) J from Brekhovskikh [Ref. 8:pp. 76-77] it can be 


shown that 


= 2Akr) Jost) KYVo) O-¥o) gy 


eee S77) 


(4.30) 


Our two solutions, given by Eqs. (4.29) and (4.30), are now in a comparable form and 


we can determine the value of the constant A(k,) to be 


ik vik 
Atk, = eeEy G0) (4.31) 


81 ky (Yo) 


Note the fact that the two terms involving the wave number in the vertical or depth 


direction, ky (y,), in Eq. (4.31) are not partially cancelled. This is due to the fact that 
their value may be either purely real (propagating wave) or purely imaginary (evanescent 
wave) caused by the infinite range of k,. As a consequence, Eq. (4.31) is the most 
concise way of expressing this quantity. 

The fact that we were able to get our solution from the WKB approximation, Eq. 
(4.30), into the same form as the Green's function, Eq. (4.29), for the free-space,isospeed 


problem suggests that our form of the WKB solution is correct. Combining Egs. (4.31), 


(4.21), and changing the variable of integration from k, to f, it can be shown that 


n= [~ eee, (k spall mie Gover (4,32) 
o ky Wo) ky @)| Sea ie aia 
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which is our solution based on the WKB approximation for a sound-speed profile that is a 
function of depth only. This is one of the fundamental results of this chapter. Now that 
this quantity is known, we can calculate the ocean medium transfer function based on the 
WKB approximation. 

The simplest way to calculate the ocean medium transfer function is to compare Eq. 
(4.32) with the overall system complex frequency response of Eq. (3.14) calculated for a 


single point source of unit amplitude. Under these circumstances MT, NT, and c,,,(f) 


equal unity and Eq. (3.14) reduces to 


Hf,r) = 2n f Hy (ffps¥q3y) Jq(2nfpt) CY IO “Yn ae (4.33) 
0 


The reason that we are interested in Eq. (4.33) is that these are the very same conditions 
that led to the solution given by Eq. (4.32). Consequently, Eqs. (4.33) and (4.32) should 
be the same. Any differences between these two equations must therefore be accounted 
for by the ocean medium transfer function as this is the only unknown quantity in the 


equations. It can therefore be shown that 


fae 
WkyVo)] 3 f kKy@ db j2ntyoX¥-Yn) (4.34) 


(a 
2ky (Vo) kyo] 


Hy Ge ae >y) = 


where ky (y) is as defined by Eqs. (4.20a) and (4.20b). Equation (4.34), then, is our 
ocean medium transfer function based on the WKB approximation and is the fundamental 


result of this chapter. 
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V. THE SIGNAL GENERATOR 


A. THE COMPLEX ENVELOPE 


Consider a transmitted electrical signal of the following form : 


x(t) =a(t)cos( 2nf,t + O(t) ) (5.1) 


This is an amplitude and angle modulated carrier where we will assume that a(t) and O(t) 
are arbitrary, real baseband amplitude and angle modulating signals, respectively. This 
implies that both a(t) and @(t) are bandlimited to the region -W <f<W. These signal 
components, a(t) and @(t), are commonly known as the natural envelope and phase of the 
signal x(t), respectively. As a consequence, x(t) is a real bandpass or bandlimited 
signal. This means that the amplitude spectrum of the signal is zero outside a frequency 
band of width 2W Hz centered about +/- f¢ Hz, where fc > W Hz and fc is referred to as 
the center or carrier frequency. 

Using Eq. (5.1) for the model of x(t), we will make use of the concept of representing 
a real bandpass signal in terms of its complex envelope [Ref. 3:pp.176-188] [Ref. 10:pp. 
74-90]. The complex envelope of x(t), denoted by x(t), is defined as 


X(t) = [ x(t) + fx(t) Jeet (5.2) 


where x(t) is the Hilbert transform of x(t). If x(t) is real, then its Fourier transform 


exhibits Hermitian symmetry. Consequently, the frequency spectrum along the negative 


frequency axis is completely determined by the spectrum for positive frequencies (and visa- 
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versa). Thus, if we remove the frequency spectrum along the negative (or positive) 
frequency axis, we will not loose any of the information contained in the signal. 


The Hilbert transform of x(t) is defined by 
x(t) = J x0 dt (5.3) 


which is the convolution of x(t) with the time function l/nt. This has the effect of 
producing a -/+ 90° phase shift for all positive/negative frequencies of the input signal. 
The Hilbert transform applies to any signal that is Fourier transformable, and by looking at 


Eq. (5.3) in the frequency domain, we obtain the analogous definition 
x(t) = IFT { -jsgn(f) X(O) } (5.4) 


where sgn(f) is the signum function and IFT{-} is the inverse Fourier transform. The 
effect of adding a real signal to its Hilbert transform is to remove the negative frequency 
components, producing a one-sided frequency spectrum which is known as the 
pre-envelope x+ (t). [Ref. 10:pp. 74-75] [Ref.11:pp. 67-72] 

This one-sided spectrum is then translated down in frequency so as to be centered at 0 
Hz to obtain the complex envelope. Thus, the complex envelope of a real bandpass signal 
is typically a complex signal with a lowpass or baseband frequency spectrum that has the 
same information content as the original bandpass signal. 

Applying the definition of the Hilbert transform to the specific case of the real bandpass 
signal x(t) in the form of Eq. (5.1) results in 
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x(t) =Xx,(t)sin(2nf-t) + x.(t)cos(2nfet) (5.5a) 

where 
X,(t) = a(t)cosO(t), (5.5b) 
x, (t) = a(t)sin@(t), (5.5c) 


and X¢ (t) and xg (t) are known as the in-phase and quadrature components of the bandpass 
signal, respectively. It should be reiterated here that one of the underlying assumptions in 
this treatment of Eq. (5.1) is that x(t) is a bandpass signal. This in turn requires X¢ (t) and 
Xs (t) to be bandlimited to the interval -W <f<W where fe > W. Equation (5.5a) is 
sometimes referred to as the canonical form of x(t). {Ref. 3:p.183] [Ref. 10:p. 85] 

By combining Eqs. (5.1), (5.2) and (5.5), it can be shown that 


HO Ow (5.6) 


is the complex envelope of the assumed form of the transmitted electrical signal given by 
Eq. (5.1). Given our initial assumptions, this means that the amplitude spectrum of the 
complex envelope described by Eq. (5.6) will be baseband centered around 0 Hz with a 
bandwidth of W Hz. 

The form of Eq. (5.6) demonstrates the advantages of complex envelope notation. 
First, it allows for the simple representation of amplitude and angle modulation. Second, 
for our assumed form of x(t) given by Eq. (5.1), which was a modulated carrier, the 


complex envelope is independent of the carrier frequency. The significance of this 
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independence from the carrier becomes apparent when dealing with the discrete 
representation of x(t), required for signal processing by digital computers. Since the 
information content of the signal x(t) is completely represented by its complex envelope 
x( t), we can sample the baseband complex envelope rather than the original bandpass 
signal. As the maximum frequency component present in x(t) is W Hz rather than 
(f¢ + W) Hz as in x(t), a lower sampling rate is required to satisfy Nyquist's criterion. 
Consequently, less tme samples or data points will be required to represent the signal 
resulting in a reduction of computation time. 

It can be shown that both the onginal bandpass signal x(t) and its Fourier transform 


X(f) can be expressed in terms of the complex envelope using 
as j2nf 
x(t) =Re[x(t)ey ©] (5.7) 
and 


X() = 5 [X(Ff) +X (-(Ff))1 (5.8) 


From the previous equations it can be seen that it is a straight forward matter to move 
between the complex envelope and common signal representations of x(t) in both the time 
and frequency domains. Thus the idea of working with the complex envelope is not at the 
expense of any substantial computational penalty in going from one representation of the 


signal to the other. [Ref. 10:pp. 85-86] 
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B. THE CW AND LFM PULSE 
For the purposes of this thesis, two classic waveforms used in underwater acoustics 
will be modelled. They are the CW (Continuous Wave) and LFM (Linear Frequency 
Modulated) pulses, both of which are real bandpass signals that can easily be expressed in 
the general form of Eq. (5.1). 
A CW waveform is simply one in which there is no angle modulation. Thus, the 
expressions for the bandpass CW signal and its complex envelope can easily be found by 


setting @(t) equal to zero in Eqs. (5.1) and (5.6) yielding 
X ow (t) = a(t)cos(2nf,t), (5.9) 
and 
Xow (t) = a(t). (5.10) 


Frequency Modulation is one of the most commonly used types of angle modulation 
and is where the instantaneous frequency, fj (t), is varied linearly with some arbitrary 
baseband signal. The instantaneous frequency is related to the time derivative of the 
instantaneous phase and, for the real bandpass signal of the form of Eq. (5.1), it can be 


shown that 


fi=fo +e e080) (5.11) 


where fc is the carrier frequency. [Ref. 10:pp. 180-183] 
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If the instantaneous frequency is a linear function of time rather than some arbitrary 
baseband signal, then Eq. (5.1) is referred to as a LFM waveform. Let the angle 
modulating signal be of the form 


@()=Dyt (5.12) 


where Dp is referred to as the phase-deviation constant or phase sensitivity. Substituting 


Eq. (5.12) into Eq. (5.11) results in the following expression for the instantaneous 


frequency : 


Dot 
fe We loess ‘ (5.13) 


The value of the instantaneous frequency given in Eq. (5.13) is a linear function of time as 
desired, composed of a constant carner frequency and the term (Dpt / ™) which is known as 
the frequency deviation. The value (Dp /7) in this instance is sometimes referred to as the 
frequency sensitivity of the modulator. [Ref. 3:pp. 193-200] [Ref. 12:pp. 267-271] 

Upon substituting Eq. (5.12) into Eqs. (5.1) and (5.6), we obtain the following 


expressions for the LFM waveform and its complex envelope : 
ifm (0) = alt)cos(2nft + Dpt’) (5.14) 


and 


~ jD ‘S 
Xum(t) = alte? Po. (5.15) 
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On comparing the above expressions with Eqs. (5.9) and (5.10), it can be seen that the 
expressions for the CW waveform can be easily obtained from those for the LFM 
waveform by setting the phase deviation constant equal to zero. This has important . 
consequences for the computer code that generates these waveforms. It means that we 
only require code for the LFM waveform since we have seen that CW is merely a special 
case of LFM. 

For the purposes of this thesis we will be concerned primarily with the rectangular 


amplitude modulating function given by 


A, |t/<$Tp 


a(t) = ie , [tl> Tp 


(5.16) 


where A is a constant and Tp is the pulse length (in seconds). It should be noted from Eq. 
(5.16) that the pulse is assumed to be centered at t = 0, and that the modulating function is 
in essence acting like a weighted on/off switch, either scaling the waveform by a constant 
or setting it equal to zero. The reason that we are interested in this particular amplitude 
modulating function is that it is because it corresponds to that transmitted by a "typical" 
SONAR system. Other amplitude modulating functions that will be considered are the 
Hanning and Hamming weighting functions, which will be discussed in the final section of 
this chapter. 

For a LFM pulse with a pulse length of Tp (in seconds), the quantity (DpTp/%) is 
referred to as the swept bandwidth (SBW) of the signal x(t). This value represents the 
frequency change (in Hz) exhibited by the signal over the period of a single pulse and 
assumes that only one sweep occurs during the period of the pulse. We are further 
assuming that this swept bandwidth is centered around the carrier frequency. 


[Ref. 3:pp. 197] 
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C. FOURIER SERIES REPRESENTATION 

AS we are modelling the ocean medium in terms of transfer functions, it is easier and 
more time efficient to work in the frequency domain as opposed to the time domain. 
Accordingly, we require x(t) to be in a form for which the transform X(f) can be easily 
computed. Except for a few simple cases, there is no closed form expression for X(f). In 
fact, X(p for the rectangular envelope LFM pulse is similar in form to a Fresnel diffraction 
integral. 


Assume that we can expand our complex envelope into the truncated complex Fourier 


series 
x(t) = S Ge. tO ee (5.17) 
q=-K 
where 
aie, * xe MO! at (5.18) 


is the complex Fourier series coefficient at harmonic q, fo = (1/To), and To is sometimes 
referred to as the data record length and is the time interval over which we are attempting to 
approximate the signal. Since this series is expanded about a period of To, the 
fundamental frequency of the Fourier series expansion is fy. Accordingly, the coefficients 
represent the magnitude and phase of the frequency components at multiples of fo Hz. 
This can be clearly demonstrated by taking the Fourier transform of Eq. (5.17) which 


yields the result: 
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K 
X() = DY) cgS(f-afy) - (5.19) 
q=-K 


The value K determines the number of terms contained in the truncated Fourier series and 
can be interpreted as the number of harmonics required to represent the complex envelope 
x(t). 

The effectiveness of the series representation given by Eq. (5.17) lies in the ability of a 
small number of terms to provide an adequate approximation of the original signal. The 
complex exponentials 


j2mq f 
eemee (5.20) 


used in this series expansion form an orthogonal basis set of functions [Ref. 15:pp. 214- 
215]. Consequently, the approximation from Eq. (5.17) using the coefficients calculated 
from Eq. (5.18) results in the minimum mean squared error estimate of x(t) for the number 
of terms used [Ref. 11:pp. 32-36]. The choice of complex exponentials as the basis set of 
functions produces a series expansion that relates the time and frequency domains. We 
have to consider a truncated series because we can only process a finite amount of data, and 
the less terms we need to represent x(t) the less computer processing will be required. 

The preceding discussion of the Fourier series assumed that x(t) was a continuous time 
signal but the concepts and results can equally be applied to a discrete time version. We 


will now discretize x(t) and introduce the Discrete Fourier Transform ( DFT{-}). The 


forward and inverse DFT can be defined as 
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? 


L : 
DFT{x@}= ¥ xe" (5.21) 


1=-L’ 
and 
‘ya j2nql/L 
DFT{Xq}=; D Xie, (5.22) 
q=K 


respectively, where IDFT{-} is the Inverse Discrete Fourier Transform and L = (2L' + 1) 
is the total number of time samples taken of x(t) during the time period Tg. To digitize 
x(t), we will sample it at time intervals of (To/L) which results in x(1) = x(ITo/L). 
Substituting this digitized version of x(t) into Eqs. (5.17) through (5.19) and comparing 
with Eqs. (5.21) and (5.22), it can be shown that 


cq => X(q) = - DFT{x(1)} (5.23) 


and 


x(l) = Lx IDFT{cq}. (5.24) 


Equations (5.23) and (5.24) in conjunction with Eqs (5.7) and (5.8) are the primary 
equations that we will use in moving between the time and frequency domains of both the 
bandpass signal and its complex envelope. 

It should be noted that this concept of representing a signal by a truncated Fourier series 
is completely general and can be applied to any signal. No assumptions regarding the 


form of the complex envelope, x(t), were made in the development of this section. 
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Consequently, if the complex envelope is of the form given by Eq. (5.6), no assumptions 


have been made regarding the amplitude and angle modulation terms of Eq. (5.1). 


D. IMPLEMENTATION 

Let us now consider the specific case of rectangular amplitude modulated CW and LFM 
pulses. We require a method of determining the value of K, so that the Fourier series 
representation contains an adequate number of harmonics to accurately represent x(t), and 
a value of L such that the sampling interval is sufficient to satisfy Nyquist's criterion or 
some other specified sampling rate. 

The complex envelope of a CW pulse can be thought of in the time domain as a 
rectangular window convolved with a sequence of delta functions at intervals of To 
seconds. It can easily be seen that as a consequence, the resulting DFT will be a discrete 
representation of the sinc waveform. To represent this signal we will take all the 
frequency components (separated by 1/To Hz) occurring within three sinc function zero 
crossings (each separated by 1/Tp Hz) from the origin. All the frequency components 
extending past this point will be less than 10% of the value of the frequency component at 
the origin, so we will consider them as being insignificant and ignore them, which yields 


the relationship 


(5.25) 


As stated earlier, there is no closed form expression for the LFM pulse in the frequency 
domain. However, for a rectangular envelope LFM pulse Papoulis [Ref. 12:pp. 267-271] 


shows that if the following condition is satisfied 
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(2f, - SBW) nae (5.26) 


then the resulting transform can be approximated as having a constant magnitude within a 
frequency band equal to the swept bandwidth, centered about the carrier frequency fc . 
Outside of this region the transform decays to zero. We need to approximate this rate of 
decay so as to decide how many frequency components we need to accurately represent the 
original bandpass signal x(t). 

The corresponding complex envelope of the LFM pulse will therefore be approximately 
constant within the region -SBW/2 < f < SBW/2 and decaying to zero outside of this 
region. Since the frequency components are spaced 1/To Hz apart, the value of K required 


to represent this constant region of the spectrum is given by 


x ~ SBWTo 
=. 


(5.27) 
Now consider the CW pulse from the point of view as being a special case of the LFM 
pulse where the swept bandwidth equals zero. We would therefore expect that in the limit 
as the swept bandwidth goes to zero that Kjfm = Kew. We will therefore approximate 
this decaying region of the spectrum as being the width of three zero crossings of the sinc 
function resulting from just the rectangular amplitude function. It is surmised that 
frequency components outside this region will be less than 10% of the frequency 
component at the origin and can therefore be disregarded. Therefore, adding the K values 


from Eqs. (5.25) and (5.27) yields 


SBWTy 3To, 
+ 
2 a 


(5.28) 








Pp 
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which is valid for both CW and LFM pulses. It should be noted from Eq. (5.17) that K is 
an integer value, so in our computer implementation of Eq. (5.28), K is rounded up to the 
next highest integer so as not to increase the truncation error. 

Like the value for K, L is also an integer, but with the additional constraint that it is an 
odd number. This is in order to satisfy the assumed form of Eq. (5.21) and to make the 
simulation output conform to that required by the adaptive beamforming algorithm that will 
be used to validate some of our results. The number of samples required can easily be 


computed from the calculated values of K using 


L=SRxK (5.29) 


where SR is the user defined value of the sampling rate and K is the value determined 
from Eq. (5.28). In order to satisfy Nyquist's criterion, the sampling rate must be greater 
thantwo. This means that the sampling frequency must be greater than twice the highest 
frequency component present in the truncated series representation of the transmitted 
electrical signal. In the computer implementation of Eq. (5.29), L is rounded to the next 
highest odd integer value in order to satisfy all of the previously mentioned constraints. 

In the computer code, the amplitude and angle modulated signal of Eq. (5.1) is defined 


by the following variables whose values are supplied by the user: 


¢ AMP : the magnitude of the rectangular amplitude modulating signal. 
¢ TP: the pulse length of the transmitted pulse (sec). 


¢ PRF: the pulse repetition frequency (Hz). This value is the same as fo and is 
equal to the inverse of the data record length To (sec). 


¢ FC: the carrier frequency (Hz). 


¢ SBW : the swept bandwidth (Hz). To describe a CW pulse, set this 
input parameter to zero. 
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* CHIRP : indicates whether the LFM sweep is increasing in frequency (up) or 
decreasing (down). 


¢ SR: sampling rate. This value is the ratio of the sampling frequency to the 
maximum frequency component present in the signal and must be greater 
than two in order to satisfy Nyquist's criterion. 


* AMOD : the amplitude modulation function, selected from a menu of choices. The 
rectangular amplitude modulation function is denoted by AMOD = 0. 


E. THE RECTANGULAR AMPLITUDE MODULATED PULSE 
To validate the theory of the preceding section, examples of rectangular amplitude 
modulated CW and LFM pulses were reconstructed from the calculated Fourier series 
coefficients of its complex envelope, denoted xce(n), and compared with the actual 
sampled bandpass signal denoted x(n). In the following plots, x(n) is plotted as a dashed 
curve while xce(n) is plotted as a solid curve. For plotting purposes, a third-degree 
polynomial curve fitting routine is used to interpolate between the data points. 
1. CW Pulse with 50% Duty Cycle 

Duty cycle is defined as the ratio Tp/To, which for this case means that the 
amplitude weighting function is nonzero for 50% of the time. While this is not a 
particularly practical transmit signal for underwater acoustic applications, it will 
demonstrate that a real bandpass signal can be reconstructed from its complex envelope. 
Combining this high value of duty cycle with a suitably low carrier frequency, the 
comparative plots of xce(n) and x(n) will be uncluttered and will help illustrate some 
important points. 

a. Pulse Reconstruction 


A CW pulse with the following characteristics is shown in Fig. 5.1a 


* AMP =1, TP =0.0125 sec, PRF = 40 Hz, FC = 400 Hz, SR = 4.0 
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It can be seen from Fig. 5.1a that the two curves, x(n) and xce(n), agree 
fairly closely. This indicates that we are indeed able to reconstruct a real bandpass signal 
from its complex envelope. However, there are some important effects to be observed. 
Note the manner in which xce(n) overshoots x(n) when the pulse is turned on and off. It 
can also be seen that Xce (n) tends to oscillate around the true value of x(n). 

A truncated Fourier series, by its very nature, will exhibit high frequency 
oscillations about the true value of the function. While these oscillations usually have 
sufficiently small amplitudes to be of no major consequence, this is not the case at a 
discontinuity. Where x(n) is turned on/off corresponds to just such a discontinuity in the 
rectangular amplitude weighting function. As more terms are added to the truncated series 
the overshoot does not tend to disappear. Instead, it becomes narrower due to the higher 
frequency terms present and the first overshoot approaches a value of 0.08949 of the 
discontinuity. This behavior is known as "Gibbs Phenomenon" and is due to the truncated 
Fourier series representation of the complex envelope, rather than the fact that we have 
used the complex envelope to reconstruct the original bandpass signal. 
{Ref. 13:pp. 107-112] [Ref. 14:pp. 179-185] 

b. Lanczos Smoothing 

It was suggested by Cornelius Lanczos [Ref. 15:pp. 219-221] that this 
phenomenon could be reduced if the Fourier series could be made to converge more 
quickly. He showed that if the series coefficients were modified by a correction term he 


called the sigma factors, given by 


sin(q7/K) 


o(K,q) = “(qniK) > (5.30) 
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where K is the order of the last term in the series, then the tendency of the high order terms 
to make a series divergent could be counteracted. As can be seen from Eq. (5.30) this was 
achieved by increasing the amount of attenuation with increasing order of the Fourier 
coefficients. Hamming [Ref. 13:pp. 107-112] comments that K in Eq. (5.30) is usually 
replaced with (K+1), and this is the convention we will adopt in our application of the 
sigma factors for this thesis. 

Lanczos also showed that the application of the sigma factors was 
analogous to averaging over an incremental time period of tT ,/2K. This means that each 
original point in the truncated Fourier series, x(t), is replaced by the arithmetic mean around 
that point, between the limits of the incremental time period. Since the ripple in the 
truncated series representation has the period of either the first term neglected or the last 
term kept, the main effects of this ripple can be removed by integrating or averaging over 
this period. The averaging interval usually chosen for this smoothing is To /(K+1), which 
corresponds to replacing K with (K+1) in Eq. (5.30). If K is large, then the region of 
local averaging is very small and has very little effect, except in the neighborhood of the 
discontinuity. This is a desirable result since the discontinuity is the region where the 
smoothing is most needed. 

When applied to the truncated Fourier series, this is known as "Lanczos 
Smoothing". Applying this smoothing to the truncated series of Eq. (5.17) yields the 
following: 
K af 


= jaraqf ot 
ae DY Odee 9 ys (5.31) 
q=-K 
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As can be seen from Eq. (5.31), the smoothed series is nothing more than the original 
Fourier series with its coefficients multiplied by the corresponding sigma factors. This 
result is what we would expect owing to the linearity of the Fourier series. The effect of 
time averaging is the same as convolving the time sequence with a weighted rectangular 
window. This in turn is analogous to multiplying the corresponding Fourier series 
coefficients of the time sequence with those of the window. Thus, the sigma factors can 
be thought of as the Fourier series coefficients of this weighted rectangular window. 
c. The Smoothed Pulse 

The effect of the application of these sigma factors can be seen in Fig. 5.1b, 
where the Lanczos Smoothing has been applied to the Fourier Series coefficients of our 
CW pulse. As expected, the oscillations and overshoot at the discontinuity have been 
reduced but at the expense of slightly broadening the pulse and reducing the rate of change 
at the leading and trailing edges of the pulse. It should be noticed that the oscillations have 
been reduced for both portions of the duty cycle, i.e., not just during the pulse. 

The first effect is to be expected since we are effectively convolving our 
time sequence with a windowing function. Whenever we convolve two sequences the 
result is a sequence larger than, but related to, the size of the original functions. Since the 
windowing function is significantly smaller than the pulse duration, we expect the effects 
of this spreading to be minimal. On comparing Figs. 5.1a and 5.1b we can see that this is 
certainly the case for this example. 

The second effect is due to the fact that the smoothing technique attenuates 
the high order Fourier coefficients. Attenuating high order coefficients is analogous to 
attenuating the high frequency components of the series representation of the bandpass 
signal. Thus, we would expect the smoothed series representation of the signal, xce(n), to 


be slower than the original Fourier series at reacting to sudden change, i.e., a discontinuity. 
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Fig. 5.la CW pulse, 50% duty cycle, no smoothing. 
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TIME T(MSEC) 


CASE: FIG 5.1B WAVEFORM:CW PULSE ( NFREQ?: 13 ) 
A: 1.0 JP: 12.5 MSEC PRF: 40.0 HZ FC: 400.0 HZ 


Fig. 5.1b CW pulse, 50% duty cycle, Lanczos smoothing. 
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Regarding the application of the sigma factors to the specific problem of Gibbs 
phenomenon, Lanczos [Ref. 15:pp. 225-226] notes that unlike the often-used Fejér 
arithmetic mean method, the sigma factors do not completely eliminate the oscillations due . 
to the discontinuity. However, the oscillations are significantly reduced and the sigma 
factors do not cause nearly as drastic a reduction in the rate of change at the leading and 
trailing edges of the pulse. It should be noted that this smoothing invalidates the 
minimum-mean-square-error property of the Fourier series. The original Fourier series 


coefficients are the only ones that satisfy this property and now they have been changed. 
2. LFM Pulse with 50% Duty Cycle 


A LFM pulse with the following characteristics is shown in Fig. 5.2a 


- AMP = 1, TP =0.0125 sec, PRF = 40 Hz, FC = 400 Hz, SRATE = 4.0 
SBW = 60Hz, "up" chirp 


The purpose of Figs 5.2a through 5.2c is to demonstrate that our analysis of the LFM pulse 
was correct. As can be seen from Fig. 5.2a, we are able to reconstruct the LFM 
waveform from the Fourier series representation of its complex envelope, as was the case 
for the CW pulse. The LFM pulse in this case is said to have “up chirp” because its phase 
deviation constant Dp is positive, that is, the frequency is increasing during the period of 


the pulse. 


As for the CW pulse, the reconstructed signal xce(n) oscillates about its true 
value x(n) and exhibits the characteristic overshoots at the leading and trailing edges of the 
pulse. The application of the sigma factors in Fig. 5.2b shows that the Lanczos smoothing 


is as effective for the LFM case as it was for the CW pulse and exhibits the same results. 
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Fig. 5.2a LFM pulse, 50% duty cycle, "up chirp", no smoothing. 
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CASE: FIG 5.28 WAVEFORN:LFM PULSE ( NPREO:15 ) 
A: 1.0 TP: 12.5 MSEC PRF: 40.0 HZ FC: 400.0 HZ SWPIBW: 60.0 HZ CHIRP:U 


Fig. 5.2b LFM pulse, 50% duty cycle, "up chirp", Lanczos smoothing. 
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It sh6uld be noted however, that since the swept bandwidth is low and the pulse 
repetition frequency is fairly high, the value of K given by Eq.(5.27) is less than one. 
This effectively means that the signal is approaching the limit of becoming a CW pulse. 
Inspection of Eq. (5.28) shows that the dominating effect on the number of significant 
frequency components is due to the rectangular pulsed nature of the waveform rather than 
the frequency modulation. Further investigation needs to be carried out on whether or not 
Eq. (5.27) is a valid saucer This will be done by the study of a LFM pulse witha 


smaller duty cycle to follow. 


The LFM pulse in Fig. 5.2c differs from the previous LFM pulses in that it is a 
“down chirp" pulse, that is its phase deviation constant is negative and the frequency is 
decreasing during the period of the pulse. The purpose of this figure is to show that the 


computer code can model the "down chirp” pulse as well as the "up chirp” case. 
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CASE: FIG 5.20 WAVEFORH:LFM PULSE ( NFREQ: 15 } 

A: 1.0 TP: 12.5 MSEC PRF: 40.0 HZ FO: 400.0 HZ SWPIBH: 60.0 HZ CHIRP:D 


Fig. 5.2¢c LFM pulse, 50% duty cycle, “down chirp", Lanezos smoothing. 
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3. CW PULSE WITH 8% DUTY CYCLE 


A pulse waveform with a small duty cycle is of more practical use for 
underwater acoustics applications, so we need to check that our computer code will handle 
more realistic signals. A CW pulse with the following characteristics is shown in Figs. 


§.3a and 5.3b: 


¢ AMP = 1, TP = 0.04 sec, PRF = 2.0 Hz, FC = 400 Hz, SR = 4.0 


As can be seen from the above parameters, the ratio Tp/To yields the value 0.08 which 
corresponds to a duty cycle of 8%. In Figs. 5.3a and 5.3b, only the curve reconstructed 
from its Fourier series coefficients, xce(n), is shown so that the plots do not appear too 


cluttered. 


From Fig. 5.3a we can observe that as expected, the unmodified Fourier 
coefficients produce a waveform that oscillates about the true value of x(n). The ringing, 
or oscillation, at the leading and trailing edges of the pulse is quite pronounced and give the 


appearance of the rectangular pulse shape spreading. 


The application of Lanczos smoothing to this case, shown in Fig. 5.3b, 
demonstrates that the oscillations about the true value have been reduced but not completely 
removed. This can be seen to be at the expense of a slight spreading of the pulse and a 
decrease in the rate of change at the leading and trailing edges of the pulse. As expected, 
these are the same effects as were observed for the CW pulse with the 50% duty cycle 


except that the extent of the ringing present in this case is more dramatic. 
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‘Fig. 5.3a CW pulse, 8% duty cycle, no smoothing. 
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Fig. 5.3b CW pulse, 8% duty cycle, Lanczos smoothing. 
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Why are we interested in this amount of ringing in the pulse representation? 
One of the outputs from this computer simulation will be the time-domain output electrical 
signal at each element in the receive array. We wish to use this information to illustrate the 
pulse distortion due to dispersion in the ocean waveguide. Consequently, we wish to 
"clean up" the transmitted pulse as much as possible so that we do not confuse distortion 
due to the Fourier series representation of the transmitted signal with distortion due to 


dispersion in the ocean waveguide. 


By comparing Figs. 5.3a and 5.3b, it can be seen that the application of 
Lanczos smoothing does produce a much "cleaner" pulse shape. It is our contention that it 
will be easier to observe the effects of dispersion on this smoothed pulse that on the 


unmodified pulse. 
4. LFM Pulse With 8% Duty Cycle 


A LFM pulse with the following characteristics is shown in Figs. 5.4a and 


5.4b: 


« AMP = 1, TP = 0.04 sec, PRF = 2.0 Hz, FC = 400 Hz, SRATE = 4.0 
SBW = 120 Hz, “up chirp" 


As for the CW case, we have only presented the signal reconstructed from its Fourier series 


coefficients, Xce(n), so that the plots do not become too cluttered. 


It can be observed that all the comments made regarding the CW pulse with the 
8% duty cycle are applicable here. It can be seen from comparing Figs. 5.3a and 5.4a 


however, that the oscillations about the true value and the amount of ringing appears to be 
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Fig. 5.4a LFM pulse, 8% duty cycle, "up chirp", no smoothing. 
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Fig. 5.4b LFM pulse, 8% duty cycle, "up chirp", Lanczos smoothing. 
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less for the LFM pulse than it does for the CW pulse. This is probably due to the number 


of frequency components used to represent the signal. 


It would appear that the value of K given by Eq. (5.28) is a more conservative 
estimate of the number of frequency components required to represent the LFM pulse than 
the estimation for the CW pulse given by Eq. (5.25). As the number of frequency 
components used is increased, the amount of ringing at the leading and trailing edges 
should be reduced. Since the amount of ringing for the LFM waveform is less than for the 
CW case, we assume that we erred on the side of caution and included more frequency 


components that absolutely necessary to represent a recognizable LFM pulse shape. This 
is consistent with the logic we used in proposing the form of Eq. (5.28) 


This would also explain why the application of the Lanczos smoothing in Fig. 
5.4b appears to have distorted the LFM pulse shape slightly less than it did for the CW 
pulse in Fig. 5.3b. If we are at the minimum number of frequency components required to 
represent the signal, then the attenuation of the high order terms will have a more dramatic 
effect than if we have a slight excess in the number of frequency components required. 
This effect can be shown quite dramatically in Figs. 5.4c and 5.4d. 

In Figs. 5.4c and 5.4d, we have reconstructed the LFM pulse using the number 
of frequency components given by Eq. (5.27), which assumes that only the constant 
portion of the spectrum of the LFM pulse is significant. As can be seen from Fig. 5.4c, 
the reconstructed pulse xce(n) is far from being a clean looking pulse shape. It oscillates 
erratically about the true value and the amount of ringing exhibited is quite extensive. 
When we apply the Lanczos smoothing to this case the waveform appears distorted as 
expected. In fact, it is distorted to the extent that we are not able to tell by looking at it that 


it has a rectangular pulse shape. 
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Fig. 5.4c LFM pulse, 8% duty cycle, "up chirp", no smoothing, 


insufficient frequency components. 
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Fig. 5.4d LFM pulse, 8% duty cycle, "up chirp", Se smoothing, 


insufficient frequency components. 
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Clearly, the assumption that only the constant portion of the spectrum of the 
LFM pulse is significant is not valid. Inspection of Figs. 5.4a through 5.4d justifies our 
scheme to calculate the number of frequency components required to represent a LFM pulse 
shape given by Eq. (5.28). Note that Eq. (5.28) was derived for the case of rectangular 
amplitude modulated pulses. Since the rectangular-shaped function contains a severe 
discontinuity, we expect that the number of frequency components given by Eq. (5.28) will 
coincide to a worst case situation. Consequently, we would expect that if we chose a 
modulating function that did not have as large a discontinuity at the leading and trailing 
edges of the pulse, we would not require as many frequency components to accurately 


represent the signal. 


F. ALTERNATIVE AMPLITUDE MODULATION FUNCTIONS 
As discussed earlier in this chapter, although we are primarily interested in the 
rectangular amplitude modulating function, it is not the only one that will be considered. 


Of interest is the Hanning weighting function given by 


: ) <1 
oe i’ [0.5 + O.Scos(2nt/Tp)] , |t]S Tp/2 (5.32) 
0 » |t]> Tp/2 
and the Hamming weighting function given by 
— p A[0.54+ 0.46cos(2nt/T,)] , |t|S Tp/2 
a(t) = { 0 salt Eeabahe? oe) 


where A is a constant and Tp is the pulse length (in seconds). It should be noted from 
Eqs. (5.32) and (5.33) that the pulse is assumed to be centered at t = 0, as was the case for 


the rectangular amplitude modulating function given by Eq. (5.16). Note that Eqs. (5.32) 
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and (5.33) are the actual pulse amplitude modulation functions and should not be confused 
with the concept of windowing functions applied to a block of data, i.e., the modulation 
functions are applied specifically to the pulse (of length Tp seconds) and not to the entire 
data record (of length To seconds). 

Inspection of Eqs. (5.32) and (5.33) reveals that they are both a maximum at the center 
of the pulse, and decrease towards zero at the leading and trailing edge of the pulse in a 
"smooth" fashion. It can also be readily seen that these two equations are very nearly 
identical. Why then are we interested in two very similar amplitude modulation functions 
of this form? From the discussion in the previous section, we expect that it will require 
less frequency components to accurately represent the transmitted signal than was the case 
when a rectangular amplitude modulation function was used. To test this hypothesis, 


examples of amplitude modulated CW pulses are plotted as follows: 


« Fig. 5.5a. Hanning amplitude modulation, 50% duty cycle, 
¢ Fig. 5.5b. Hanning amplitude modulation, 8% duty cycle, 

¢ Fig. 5.6a. Hamming amplitude modulation, 50% duty cycle, 
¢ Fig. 5.6b. Hamming amplitude modulation, 8% duty cycle. 


In Figs. 5.5a and 5.6a, the pulse reconstructed from the calculated Fourier series 
coefficients of its complex envelope, denoted x¢e(n), is compared with the actual sampled 
bandpass signal, denoted x(n), where x(n) is plotted as a dashed curve while xc¢e(n) is 
plotted as a solid curve. For plotting purposes, a third-degree polynomial curve fitting 
routine is used to interpolate between the data points. In Figs. 5.5b and 5.6b, just xce (n) 


is plotted to prevent the plot from becoming too cluttered. 
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Comparing Figs. 5.5a and 5.6a with Fig. 5.la, it can be seen that indeed fewer 
frequency components are required to accurately represent the transmitted signal when the 
amplitude modulation function is of the form given by Eqs. (5.32) or (5.33), than if it is of 
the form given by Eq. (5.16). Whereas we introduced the concept of Lanczos smoothing 
for the rectangular modulated pulse in an attempt to obtain a "cleaner" representation of the 
signal, we find that this is totally unnecessary for the Hanning and Hamming modulation 
functions. The observations that can be made when comparing Figs. 5.5b and 5.6b with 
Fig. 5.3a are consistent with those just described. From these comparisons we can make 
the statement that the Hanning and Hamming modulated CW pulses require roughly half 
the number of frequency components for accurate representation than does the rectangular 
modulated pulse. No comparison of modulated LFM pulses was carried out as it was 
judged past the scope of interest for this thesis. 

Further comparison of Figs. 5.5 and 5.6 shows that for a given number of frequency 
components the representation of the Hamming amplitude modulated pulse is better than the 
Hanning modulated pulse. This means that there is closer agreement between the pulse 
reconstructed from the calculated Fourier series coefficients of its complex envelope and the 
actual sampled bandpass signal when the Hamming amplitude modulation function is used. 
In particular, the Hamming modulated pulse displays significantly less ringing before/after 
the leading/trailing edges of the pulse for the lower duty cycle case shown in Figs. 5.5b 
and 5.6b. 

The advantage in reducing the number of frequency components required to represent a 
signal is a simple one. Since the computer simulation must compute a value of the overall 
system complex frequency response for each frequency of interest, and this operation 
comprises the majority of the computational effort, a decrease in the number of required 


frequencies will mean an almost proportional drop in computation time, i.e., if the 
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Fig. 5.6b CW pulse, 8% duty cycle, Hamming amplitude modulation. 
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number of frequencies is reduced by half, then the program run time will be approximately 
reduced by half. 

However, this is not the main reason that we consider amplitude modulation functions 
of this form. As stated earlier, we are primarily interested in the rectangular amplitude 
modulating function. Our interest in the Hamming and Hanning functions stems from the 
fact that many currently available computer models dealing with sound propagation in the 
ocean deal with transmitted waveforms having these particular amplitude modulation 
functions, presumably in an attempt to reduce computation times. If we then also have the 
option of using these same modulation functions, we will be in a better position to compare 
our results and performance with other ocean medium propagation models, e.g., SAFARI 


[Ref. 7] . 
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VI. COMPUTER SIMULATION RESULTS FOR THE 
ISOSPEED MEDIUM CASE 


A. THE OVERALL SYSTEM COMPLEX FREQUENCY RESPONSE FOR AN 
ISOSPEED MEDIUM 

The first step in validating the developments carried out in Chapters 2 through 4 is to 
apply the analysis to the simplest scenario, i.e., the free-space, isospeed medium case. 
Since the preceding development of the overall system complex frequency response given 
by Eq. (3.14) assumed a free-space problem, all we need to do is apply the condition of an 
isospeed medium to the ocean medium transfer function given by Eq. (4.34). As 
previously discussed, the phase integral of the ocean medium transfer function has a closed 


form expression in the isospeed case, which results in the relationship 


SIKYU MY 0) gi2tf yo - Yn) 


2 ky Vo) ky | (6.1) 


Hy (ff. yosy) = 


where 


1/2 
PGy== om (ife,)-t] . f& S(fic.). (6.2a) 
” Wie 2: 2 
i) == i2wi (fig) ] «fr > (Phey)* (6.2b) 


ky (y,) =2nfy(y,) =2nfy, and co is the constant speed of sound in the isospeed 


medium. Equation (6.1), then, is our ocean medium transfer function for a free-space, 
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isospeed medium. Note that the ocean medium transfer function given by Eq. (6.1) 
contains no integral operations due to there being a closed form expression for the phase 
integral of Eq. (4.34). This result is significant since numerical integration routines are 
very time consuming. Therefore, we expect that the computational effort required for the 
isospeed case will form a baseline against which to compare other cases that do not have 
closed form solutions. 

It can be easily seen from Egs. (6.2a) and (6.2b) that for an isospeed medium, the 
propagation vectors in the Y directions, ky (yg) and ky (y), are equal. The ray acoustics’ 
interpretation of this is that rays travel in straight lines and will not bend. It should be 
remembered from the development in Chapter 3, that Eq. (6.2a) represents propagating 
waves while Eq. (6.2b) represents evanescent waves. 

Returning to our expression for the overall system complex frequency response given 


by 


M5 
M3 


H(f,r) = 2n f Cmn(f) Hy (of, ¥03y) Jo (20£ 5m) 
Qo m 


: 
2 


x ePMYO Ye ae (3.14) 


we notice that the region of integration with respect to f; is from zero to infinity. Equation 
(3.14) is referred to as a full-wave solution for our ocean medium propagation problem, 
and can be thought of as an infinite summation of ray acoustic solutions [Ref. 7:p. 36]. 
From Eqs. (6.2a) and (6.2b), it can be seen that this region includes both propagating and 
evanescent waves. Since the propagation vector in the Y direction has different forms 
depending on whether the wave is propagating or not, we will split Eq.(4.33) into two 


integrals. One integral will cover the region of propagating waves, i.e., integrate with 
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respect to fr from zero to (f/ Co), while the second will cover the region of evanescent 
waves, i.€., integrate with respect to fr from (f / co) to infinity. 

As discussed in Chapter 4, our general solution for the case of a unit amplitude point 
source in an unbounded, homogeneous medium should collapse into the form of the free- 
space Green's function for an isospeed medium given by 


-ikR 
-e! 


4nR 





g(r) = (4.25) 


In this chapter, we will refer to Eq. (4.25) as the exact solution for the free-space, isospeed 
medium since no simplifying approximations were made to obtain this result. Examination 
of Eq. (4.25) shows that under these circumstances, the propagating wave undergoes 
spherical spreading, i.e., the amplitude of the acoustic field is inversely proportional to 
radial distance between the transmitter and the field point. This radial distance is 
commonly referred to as the Range-Line-Of-Sight (RLOS). 

Reviewing the development leading up to the expression for the overall system complex 


frequency response due to a unit amplitude point source given by 


H(f,r) = 2n - Hu (fifesYosy) Jo(2mferm) een YOOO “Ye ae, (4.33) 
0 


remember that it is mathematically the same as the acoustic field given by Eq. (4.25). Asa 
result, the overall system complex frequency response should behave the same as the 
acoustic field due to a point source. Accordingly, if we plot RLOS versus the magnitude 
of the overall system complex frequency response on log-log axes for a particular 


frequency component of the transmitted signal, we would expect the result to be a straight 
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line. Further, due to the isospeed nature of the medium, we would expect all propagating 
frequency components to behave the same since we are only concerned in this simulation 


with the effects of spreading and have ignored losses due to absorption. 


B. RESULTS OBTAINED BY IGNORING EVANESCENT WAVES 

An initial investigation was carried out for an isospeed medium with a sound speed of 
1475 m/sec. For a transmit signal, we used the waveform of Fig. 5.1b, which was a 
Lanczos smoothed, rectangular amplitude weighted CW pulse with the following 


parameters : 


¢ AMP = 1.0, TP = 12.5 msec, PRF = 40.0 Hz, FC = 400.0 Hz 


* minimum frequency = 160.0 Hz, maximum frequency = 640 Hz. 


For this scenario, the longest wavelength (which corresponds to the lowest frequency 
component present) is equal to 9.219 m, while the shortest wavelength (which corresponds 
to the highest frequency component present) is equal to 2.305 m. Since the evanescent 
waves are exponentially decaying instead of being oscillatory, they are considered as being 
non-propagating and are usually ignored at distances greater than a couple of wavelengths 
from their source. For the case just described, we would expect evanescent waves to be 
insignificant at ranges greater than say fifty meters. Applying this assumption to our 
expression for the overall system complex frequency response, i.e., ignoring the 


integration over the region of evanescent waves, results in 


ELCs 2 
H(f,r) = 2n il Hy (Cf, ¥oy) Jo@umpe: 0° = tdi men) 
0 
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which is an approximation of the full-wave solution since the region of integration is no 
longer infinite. We will refer to any overall system complex frequency response that has 
such a limited region of integration as our approximate (full-wave) solution to the ocean 
medium propagation problem. 

It should be stressed at this point that the computer implementation of the overall 
system complex frequency response has been kept as general as possible. This means that 
instead of having to program separate forms of Eg. (3.14) for different scenarios, e.g., Eq. 
(6.3), just the generalized form of the overall system complex frequency response given by 
Eq. (3.14) is implemented in the computer code. This generalized form can then be 
modified as required, by a number of user defined variables. To obtain the overall system 


complex frequency response given by Eq. (6.3), we merely set the user defined variables 


* transmit array : MT =1, NT = 1, STEER = FALSE 
* integration routine : RATIO = 1.0 
* sound speed profile : NGRAD =0 


where MT and NT describe the dimensions (number of elements) of the transmit array, as 
discussed in eniet 2, and STEER is a logical variable that designates whether any 
beamsteering is to be done. If no beamsteering is to be done, then the transmit array 
complex weights, given by Eq. (2.25), are set equal to unity. The variable, RATIO, 
defined as 


RATIO = f,c(yo)/f , (6.4) 
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determines the upper limit of the integration with respect to fr , while the variable NGRAD 
set equal to zero designates an isospeed medium and the program uses the closed form 
expression for the phase integral appearing in Eq. (3.14). 

To demonstrate the validity of the preceding discussion, we need to determine if our 
derived overall system complex frequency response does in fact reduce to the isospeed, 
free-space, Green's function and show the expected (1 / RLOS) fall-off in the magnitude of 
the acoustic field. In other words, we need to compare the results given by the 
approximate full-wave and exact solutions. Accordingly, the magnitude of the overall 
system complex frequency response versus range-line-of-sight is plotted in Fig. 6.1 for 
three different frequencies. The frequencies chosen correspond to the minimum, carrier, 
and maximum frequencies of the transmitted signal. Additionally, the approximate full- 
wave and exact values of the magnitude of the acoustic field at the carrier frequency of 400 
Hz and for the values of RLOS used in Fig. 6.1 are detailed later in Table 6.3. 

As can be seen from Fig. 6.1, the magnitude of the acoustic field has been determined 
for values of RLOS ranging from ten meters to ten thousand meters. This is equivalent to 
a range of one to thousands of wavelengths, depending on which of the frequency 
components we are considering. From our preceding discussion, we would expect that 
the acoustic field calculated using the approximate full-wave solution, given by Eq. (6.3), 
will not reduce exactly to the isospeed, free-space, Green's function for ranges less than 
about fifty meters. This is because we have assumed that evanescent waves have a 
noticeable effect in this region, which Eq. (6.3) ignores. Conversely, would we expect 
close agreement for ranges greater than fifty meters as we have assumed that the effect of 
evanescent waves in this region is negligible. In terms of the log-log plot of magnitude 
versus RLOS shown in Fig. 6.1, we would expect the graph to be a straight line for ranges 


greater than fifty meters, and to deviate away from this ideal for shorter ranges. 


84 


Acoustic Field 
Velocity Potential ( m**2 / sec ) 


Slope of theoretical 1/r fall off 
due to spherical spreading 





RLOS (m) 


Fig. 6.1 Magnitude of the acoustic field versus RLOS when 
the effects of evanescent waves are ignored. 
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It can readily be observed from Fig. 6.1 that for long ranges, the graph does indeed 
approximate a straight line with the expected slope, and for short ranges, it deviates from 
this straight line. These results appear, at first glance, to be what we expected. However, 
note the range at which the plotted curve starts to approximate a straight line. This range is 
about five hundred meters, or ten times the range that we thought would be significant. 
Alternatively, this range represents between 54 and 216 wavelengths, depending on which 
of the three frequency components we are considering. This would suggest that the effects 
of evanescent waves are indeed still significant at ranges greater than a couple of 
wavelengths from their source. 

Consider also, the actual values of the approximate acoustic field that are contained in 
Table 6.3. Even at a range of ten thousand meters, the approximate full-wave and exact 
solutions calculated for the carrier frequency, still differ by almost 1%. For the carrier 
frequency, this range corresponds to 2,712 wavelengths, which is well beyond where we 
expected evanescent waves to have an effect. 

Close inspection of Fig. 6.1 does indicate, however, that the low frequency term, i.e., 
160 Hz, is the one that is behaving the most erratically. This is what we would expect 
with regards to the deviation from the straight line predicted by the exact solution since this 
frequency corresponds to the longest wavelength. Consequently, for any given distance, 
this is the evanescent component that would have decayed the least. Since we are ignoring 
evanescent waves in this case, this is the frequency component that should deviate most 
from the exact result. Based on these observations and trends, we will assume that our 
derivation of the overall system complex frequency response has not been proved in error 
so far, and that it is the assumption regarding the significance of evanescent waves that is 


incorrect. 
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Even though our approximate solution results deviate from those predicted by the exact 
solution, is the amount of deviation really significant ? Remember that one of the outputs 
that we wanted from the computer simulation was the time-domain electrical signal at each 
element in the receive array, which would be used to show pulse distortion due to the 
effects of dispersion, and as input data to space-time signal processing algorithms. To 
demonstrate this pulse distortion, we need to compare the original transmitted electrical 
signal with the received electrical signal and observe the changes that propagation through 
the ocean medium has caused. Such a comparison can be made using Figs. 6.2a and 
6.2b. Note that 6.2b is the received electrical signal ata RLOS of 100 m._ Referring to 
Fig. 6.1, we see that this value of RLOS is one for which the calculated acoustic field 
varies significantly with frequency. This particular RLOS was chosen to represent a near 
worst Case situation and will demonstrate the significance, if any, of the inexactness of the 
approximate full-wave solution when evanescent waves are ignored. 

It can be seen from comparing Figs. 6.2a and 6.2b, that the output electrical signal is 
indeed significantly different from the transmitted signal. Recall that there should be no 
dispersion in an isospeed, free-space medium. This means that the received electrical 
signal should be a scaled replica of the transmitted electrical signal, i.e., identical in pulse 
shape. The problem with the approximate full-wave solution that produced the output 
electrical signal shown in Fig. 6.2b is that, if we repeated the problem for a non-isospeed 
medium in which we do expect pulse distortion due to dispersion, we would not know 
how much of the distortion was due to the dispersion and how much was due to the 
approximate nature of our solution. 

Next, consider the previous problem repeated for a RLOS of 1,000 m. From Fig. 6.1 
it can be seen that for this range, the acoustic field behaves a lot more as predicted and is 


less dependent on frequency. This means that there is much closer agreement between the 
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Fig.6.2a Transmitted electrical signal. 
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Fig.6.2b Received electrical signal, RLOS = 100 m, when the 


effects of evanescent waves are ignored. 
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approximate full-wave and exact solutions than existed for RLOS = 100 m. The output 
electrical signal for this case is shown in Fig 6.2c, where the transmitted electrical signal is 
the same as for the previous case shown in Fig. 6.2b. As expected, the output electrical 
signal indeed looks a lot more like the transmitted electrical signal than was the case in Fig. 
6.2a. In fact the amount of distortion can be considered imperceptible for RLOS = 1,000 
m, and, for all practical purposes, Fig. 6.2c is a scaled version of the transmitted signal. 
Based on these results for the received time-domain electrical signal, it appears that the 
approximate full-wave solution is adequate for demonstrating the effects of pulse distortion 
for signals containing no frequency components below 160 Hiz, where the receive element 
is at a range greater than 1,000 m from the source. This situation should be acceptable for 
long range propagation problems, but for any short range calculations, it is obviously 


inadequate. 
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Fig. 6.2c Received electrical signal, RLOS = 1,000 m, when 


the effects of evanescent waves are ignored. 
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Consider next, the other use of the time-domain electrical signal at each element in the 
receive array, that is, as input data to space-time signal processing algorithms. As stated 
earlier, the computer code developed for this ocean medium propagation problem does not 
perform any signal processing on the received electrical signal at each point in the array, but 
rather, the data is written to a file so that it can be operated on by ancillary programs. Asa 
further validation of the generated results, we can operate on the time-domain data with a 
beamforming algorithm that localizes the position of the source, and compare the computed 
and actual source position. The algorithm used was an adaptive beamforming and non- 
linear estimation algorithm developed by Breckenridge [Ref. 16] that localizes a source in 
range, depression angle, and azimuthal angle, according to the geometry outlined in Fig. 
6.3. From Fig. 6.3 it should be noted that the localization of the source, calculated by this 
adaptive algorithm, is relative to the center element in the receive array. 

The data generated from the two cases, RLOS = 100 m and RLOS = 1,000 m, was 
processed with this algorithm using the output electrical signal from each of the elements in 
a 3x3 receive array, and the results are shown in Table 6.1. The values RLOS (X) and 
RLOS (Y) detailed in the table correspond to estimates of the RLOS calculated from 
information along the X and Y axes, respectively. From the results contained in Table 
6.1, we can see that the estimated values from the beamforming algorithm are closer to the 
true values for the case RLOS = 1,000 m than for the case RLOS = 100 m, which is 
consistent with our previous observations in this section. 

An interesting feature of the results is the way that RLOS (Y) differs radically from the 
true value of RLOS, while RLOS (X) is fairly close to the true value. Note that this occurs 
for both values of RLOS considered, that is, RLOS (Y) does not improve as the value of 
RLOS is increased. This phenomena is in contrast to the trends observed earlier in this 


section. The fact that the effects of evanescent waves are being ignored means that we are, 
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Fig. 6.3 Geometry for adaptive beamforming algorithm. 


RLOS RLOS (X) | RLOS(Y)| © e y wy 
(m) (m) (m) (deg) |) (deg) (deg) | Cdeg) 
true est est true est true est 
100 115.02 29.54 | 30.000 | 28.274 | 270.00 | 270.00 
1000 960.71 284.619 | 30.000 | 29.942 | 270.00 | 270.00 


Table 6.1 Adaptive beamforming results when 
evanescent waves are ignored. 
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in essence affecting the wave solution in the Y direction to a much greater degree than in the 
X direction. Recall that our form of the overall system complex frequency response, given 
by Eq. (3.14), was developed assuming a sound-speed profile that was a function of 
depth. In our development we showed that this meant Eq. (3.14) was axisymmetric. If 
the overall system complex frequency response is axisymmetric, then the relative behavior 
of the acoustic field in the X direction should be minimally affected by the behavior of the 
acoustic field in the Y direction. Consequently, if there is a problem with the description 
of the full-wave solution in the Y direction, as in this case, the value of RLOS (X) should 
still be fairly accurate, which is the result that we observed. 

These results also demonstrate that while the time-domain output electrical signal from a 
single receive array element may look correct to the eye and suggest that the accuracy of the 
generated data is adequate, it may in fact not be accurate or correct enough for signal 
processing purposes. This suggests that observing the output electrical signal is not a 
good way of validating the performance of the approximate full-wave solution. 

The remaining output that we wanted from this computer simulation was the magnitude 
and phase of the complex acoustic field incident upon a receive array as a function of 
frequency and spatial coordinates. A convenient way of viewing these results, apart from 
a three-dimensional plot, is to tabulate the magnitude and phase of the acoustic field in the 
same order as the elements of the receive array. Consider a 3x3 element receive array. 
We will adopt the element designation scheme shown in Fig. 6.4 to number the elements 
(m,n ) contained in the array. 

The values of the magnitude and phase (in degrees) of the acoustic field are shown 
in Tables 6.2a and 6.2b for the cases RLOS = 100 m and RLOS = 1,000 m, respectively. 
Note that the first value is the magnitude and the second value is the phase. Given the 


placement of the receive array relative to the source, i.e., XR =0 and YR = 1025 m, which 
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Fig. 6.4 Element designation for a 3x3 receive array. 


is below the source, we would expect the magnitude of the acoustic field to decrease with 
increasing n. This corresponds to reading down the table, or equivalently, increasing the 
depth of the receive element or field point. We would also expect the value of the acoustic 
field to decrease with increasing absolute value of m. This corresponds to reading across 
the table out from the center, or equivalently, increasing the displacement in the X direction 
of the element or field point. We expect the magnitude of the acoustic field to decrease in 
both of these situations since their effect is to increase the radial range between the source 
and the receive element , or field point. 

From these two tables it can be seen that all the values for the acoustic field do not 
behave as expected. In Table 6.2a, the first two rows of results have the magnitude of the 
acoustic field increasing as we move away from m=0. In Table 6.2b on the other hand, it 
is the third or final row of values that behaves in this fashion. In both cases however, the 


magnitude of the acoustic field decreases as the depth is increased, as expected. As for the 
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8.69822 e-4 
- 168.608 


8.47479 e-4 
139.860 


7.82107 e-4 
84.800 


8.68689 e-4 
- 167.994 


8.47246 e-4 
140.423 


7.82997 e-4 
85.376 


8.69822 e-4 
- 168.608 


8.47479 e-4 
139.860 


7.82107 e-4 
84.800 


Table 6.2a Acoustic field at array elements, RLOS = 100 m. 


8.00331 e-5 
169.065 


7.97398 e-5 
113.082 


7.93029 e-5 
56.758 


8.00375 e-5 
169.130 


7.97420 e-5 
113.150 


7.93010 e-5 
56.826 


8.00331 e-5 
169.065 


7.97398 e-5 
113.082 


7.93029 e-5 
56.758 


Table 6.2b Acoustic field at array elements, RLOS = 1,000 m. 
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results obtained from the use of the adaptive beamforming algorithm, observing the values 
of the complex acoustic field as a function of frequency and spatial coordinates alerts us to 
the fact that the data produced by our approximate full-wave solution is inaccurate. 

From the results of this section, we can formulate two primary conclusions. First, the 
assumption that the effects of evanescent waves past a distance of a couple wavelengths is 
insignificant is an invalid one, at least for the situation we are modelling and for the full- 
wave method used. Second, comparing the transmitted and received electrical signals for 
the free-space, isospeed medium case is not sufficient to determine whether or not the 
approximate full-wave solution, is accurate. To gauge the performance of the approximate 
solution we need to either observe the acoustic field as a function of spatial coordinates for 
an appreciation for whether the solution is behaving correctly, or apply signal processing 


techniques to obtain some quantitative measure for the accuracy of the solution. 


C. RESULTS OBTAINED BY INCLUDING EVANESCENT WAVES 

To include the effects of evanescent waves, we need to extend the region of integration 
with respect to fr in Eq. (6.3). As was previously discussed, the region of evanescent 
waves extends from fy equal to (f/ co) to infinity. But is it necessary to integrate over this 
entire region 9 After some trial and error, it was found that integrating with respect to f; up 
to a limit of 1.2 (f / cg) was adequate to give very good results down to a RLOS of 10 m. 
Recall that evanescent waves have more effect at short ranges. Therefore, if this upper limit 
of integration is sufficient for a RLOS of 10 m, then it will also be sufficient for longer 
ranges. Since no practical reason could be thought of as to why we would be interested in 
ranges less than 10 m, this is the upper limit of integration that we adopted and is that used 
to obtain the results of this section. Therefore, our approximate full-wave solution is now 


given by 


oe 


L2(iicy) oe 
H(f,r) = 2n i Hy (ff, .¥o:y) Jo Qmfptm) ee YOO £, df, (6.5) 
0 


Recall from our previous discussion that the region of integration in Eq. (6.5) is, for 
computational purposes, divided into two regions to account for the change in the form of 
the propagation vector in the Y direction, ky (y), when fr changes from being less than to 
greater thanf/cg. This empirically determined value for the upper limit of integration 
agrees with observations made by Schmidt [Ref. 7:p. 31] with regard to wavenumber 
integrations. 

For the analysis contained in this section, the scenario is the same as for the previous 
case where we ignored the effects of evanescent waves, i.e., same transmit signal, array 
geometry and ocean medium. To allow for this extended region of integration, all that has 
to be done is to set the variable RATIO equal to 1.2 instead of 1.0. This demonstrates one 
of the advantages in having kept our initial analysis and computer code implementation as 
general as possible. It is a very simple matter to alter the problem that requires no 
programming changes. 

The magnitude of the overall system complex frequency response was determined for a 
range of values for RLOS and for three frequencies corresponding to the minimum, carrier, 
and maximum frequency components of the transmitted signal. It was found, however, 
that, unlike the previous case, there was no discernable frequency dependence in the values 
for the magnitude of the overall complex frequency response determined from the 
approximate full-wave solution given by Eq. (6.5). As was previously discussed, this 
frequency independence is the result predicted from the exact form of the isospeed, free- 
space Green's function. Since the results for the three frequency components were 


identical, only the results for the carrier frequency, 400 Hz, are shown in the log-log plot 
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of acoustic field versus RLOS shown in Fig. 6.5. Additionally, the approximate full-wave 
and exact values of the magnitude of the acoustic field for this frequency and values of 
RLOS used in Fig. 6.5 are detailed in Table 6.3. 

From Fig. 6.5, it can be observed that the log-log plot of the acoustic field magnitude 
versus RLOS is identical in slope to the straight line predicted by the isospeed, free-space, 
Green's function, over the entire range of values of RLOS considered. Comparing Figs. 
6.1 and 6.5, it can easily be observed that there has been a radical improvement due to the 
inclusion of a small region of evanescent waves. The results from Fig. 6.5 suggest that 
the upper limit of integration chosen in the approximate full-wave solution given by Eq. 
(6.5) is adequate for values of RLOS down to at least 10 m. 

Inspection of Table 6.3 shows that with the inclusion of this region of evanescent 
waves, the magnitude of the acoustic field given by the approximate full-wave solution in 
fact agrees with the exact values predicted by the isospeed, free-space, Green's function to 
within five significant figures. The comparison was limited to this number of significant 
figures since the numerical integration was only carried out to sufficient terms to give a 
maximum of a 0.01% relative error in the value of the integral. Alternatively stated, the 
value of the acoustic field calculated from Eq. (6.5) is accurate to within +0.01%. 

Let us now consider the displaying of the received electrical signal for the purposes of 
showing pulse distortion due to dispersion. Comparing Figs. 6.6a and 6.6b, it can be 
observed that there is no discernable pulse distortion, which is what theory predicts for 
wave propagation in a free-space, isospeed medium. This result was also expected from 
experience with the analysis of the previous case, in which we ignored the effects of 
evanescent waves. In that case, even though there was some frequency dependence in the 
value of the acoustic field calculated using the approximate full-wave solution of Eq. (6.3), 


there was very little discernable distortion in the received electrical signal. For the acoustic 
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Fig. 6.5 Magnitude of the Acoustic field versus RLOS when 
the effects of evanescent waves are included. 
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RLOS (m) 1/4nR | H(f,r) | | H(f,r) | 

(no evanescent (incl. evanescent 

waves) waves ) 
10 7957 1e-3 7.9702 e-3 7.9577 e-3 
50 1.5915 e-3 1.3983 e-3 1.5915 e-3 
100 7.9577 e-4 8.4725 e-4 7.9577 e-4 
200 3.9789 e-4 3.9676 e-4 3.9789 e-4 
300 2.6526 e-4 2.5923 e-4 2.6526 e-4 
500 1.5915 e-4 1.5663 e-4 1.5915 e-4 
1,000 7.9577 e-5 7.9742 e-5 0.95 17e-5 
2,000 3.9789 e-5 3.9091 e-5 3.9789 e-5 
3,000 2.6526 e-5 2.5978 e-5 2.6526 e-5 
5,000. 1.5915 e-5 1.5753 e-5 1.5915 e-5 
10,000 7.9577 e-6 8.0182 e-6 7.9577 e-6 


Table 6.3 Comparison between exact and predicted 


magnitude values of the acoustic field. 
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field calculated using the more accurate approximate full-wave solution given by Eq. (6.5), 
we would expect our received electrical signal to look even more like a scaled version of the 
transmit signal. This was, in fact, the observed result. 

The beamforming algorithm was then applied to the output electrical signals from a 3x3 
receive array, as in the previous case, with the results being summarized in Table 6.4. The 
first thing to notice about these results is that unlike the previous case that ignored 
evanescent waves, both RLOS (X) and RLOS (Y) give accurate values for the radial 
distance between the source and the center of the receive array. Additionally, the estimated 
position angles are very accurate. 

Finally, an inspection of the complex acoustic field at the elements of the receive array 
is shown in Table 6.5. It can be clearly seen that the values of the magnitude conform to 
the expected pattern, as described for the previous section. Note that the change in 
magnitude, reading down the table, is much greater than when reading across the table. 
While this may appear to be an error at first, it is in fact correct, and is due to the 
positioning of the array. Recall that the receive array is positioned such that XR = 0, 
YR = 1050 m, i.e., the receive array is 50 m below the transmit array, and RLOS = 100 m. 
Also recall that the distance between elements in the X and Y directions is the same. 
Accordingly, it can be easily shown that the change in radial distance between the source 
and a receive element is much larger for an incremental change in the Y direction than for 
the same change in the X direction. As a result, the change in the magnitude of the 


acoustic field will be similarly affected. 


From the results of this chapter so far, we have seen that it is not a valid assumption to 
ignore the effects of evanescent waves when developing an approximation of the full-wave 
solution to an ocean medium propagation problem. However, it appears that we only have 


to include a small portion of the evanescent region for an approximate full-wave solution to 
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Fig. 6.6a Transmitted electrical signal. 
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Fig. 6.6b Reccived electrical signal, RLOS = 100 m, when the 


effects of evanescent waves are included: 
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be valid. From these results we will further assume that the region of integration with 


respect to fr , given by 


O<f,< 1.2(f/co), (6.6) 


for the overall system complex frequency response, given by Eq. (3.14), is adequate for 
sound-speed profiles other than the isospeed case investigated in this chapter. We will 
utilize this assumption in the following chapter, when we investigate the behavior of the 
acoustic field in a free-space medium characterized by a linear sound-speed profile with a 
single gradient. We have also seen that the accuracy of the approximate full-wave solution 
is a function of the numerical integration used to calculate the overall system complex 
frequency response. There is a trade-off between computation time and accuracy that will 


be addressed in the final section of this chapter. 


D. PLOTTING OF THE 3-D COMPLEX ACOUSTIC FIELD 

So far in this chapter, all our validation has taken place using a similar geometry, i.e., 
XR =Qand YR> YT. This means that the receiver has always been deeper than the 
transmitter, and has had no offset in the X direction. To have any confidence that the 
theory and computer implementation of the overall system complex frequency response is 
correct, we need results for every permutation of the X and Y direction offsets of the 
receiver with respect to the source. In other words, we need to place the receiver in a 
variety of positions, e.g., positive X offset and above the source, to fully test the 
algorithm. 

Since one of the required outputs is the magnitude and phase of the complex acoustic 
field at the locations of the array elements for use by matched-field processing algorithms, a 


3-D plotting routine was implemented and added as a visual aid to help conceptualize what 
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RLOS RLOS (X) | RLOS(Y)| © e y y 


(m) (m) (m) (deg) | (deg) | (deg) | (deg) 
true est est true est true est 
100 99.984 99.977 | 30.000 | 30.004 | 270.00 | 270.00 


Table 6.4 Adaptive beamforming results when 
evanescent waves are included. 


n/m 1 0 -1 
-] 8.00292 e-4 8.00346 e-4 8.00191 e-4 
-167.604 -166.952 -167.604 
0 7.95722 e-4 7.95774 e-4 7.95722 e-4 
136.639 137.287 136.639 
1 7.91125 e-4 7.91177 e-4 7.91125 e-4 
79.909 80.554 79.909 


Table 6.5 Acoustic field at array elements, RLOS = 100 m. 
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the acoustic field looked like. This routine plots, on separate graphs, the magnitude and 
phase of the overall system complex frequency response for a particular frequency versus 
the relative position of each of the array elements. Recall from our earlier development that 
the value of the frequency spectrum of the acoustic field at any point is just the value of the 
overall system complex frequency response times the frequency spectrum of the transmitted 
electrical signal. By then looking at the plots, particularly those of the magnitude of the 
acoustic field, it can easily be determined if the field exhibits the appropriate orientation and 
curvature for the given relative positions of the source and the receive array. 

Note that the planar array of point sources we are calling the transmitter does not 
necessarily have to be a physical transducer array. It could represent sampled values of the 
acoustic field at some spatial location (or grid of locations) due to some other source. In 
this situation, the array could be thought of as detecting an acoustic signal and reradiating 
that acoustic signal without alteration. 

As was previously discussed, we need to generate a series of plots of the complex 
acoustic field for various combinations of source / receiver locations. A series of plots 
were generated for the six different permutations of the X and Y offsets of the receive array 
relative to the source. For all these plots, the source was located at a depth of 1,000 m. 
Recall that the transmit array (or single point source in this case) is centered within the 
coordinate axis system such that XT and ZT are both equal to zero and YT is the depth 
below the sea surface (see Fig. 2.3). The receiver was a9 x9 planar array of point 
sources, which resulted in the computation of the overall system complex frequency 
response for 81 different spatial locations. These values of the overall system complex 
frequency response were calculated for a frequency of 400 Hz, which is the carrier 
frequency of the transmitted electrical signal shown in Figs. 6.2a and 6.6a. The RLOS, 


i.e., the radial distance between the source point and the center of the receive array, was 
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assumed to be only 50 m in order to reduce computer processing time. This short range 
was chosen so that we could process data a reasonable sized array which would allow us 
to observe curvature in the acoustic field, i.e., so that the plots would not just look like flat 
planes. The array size was chosen as being adequate to produce a 3-D plot, for the RLOS 
as previously chosen, in which the orientation and curvature of the acoustic field values 
was easily discernable. A constant speed of sound of 1475 m / sec was assumed for all 
cases. In the following 3-D plots, the X and Y offsets AY = YR- YT and AX = XR are 
used to describe the relative positions of the centers of the transmit and receive arrays. 

A series of six sets of plots were generated, showing both the magnitude and phase of 


the overall system complex frequency response, for the following positions of the receive 


array : 
¢ Fig.6.7 YR=1000m, XR=0m, (AY =0, AX =0) 
- Fig.6.8 YR=1050m, XR=0m, (AY =50, AX =0) 


- Fig.6.9 YR=1010m, XR=20m, (AY = 10, AX = 20) 
° Fi 


iy 


g.6.10 YR=980m, XR=10m, (AY =- 20, AX = 10) 
¢ Fig.6.11 YR =1020m, XR=- 20m, (AY = 20, AX =- 20) 
* Fi 


_ 


g.6.12 YR=990m, XR=- 20m, (AY =- 10, AX =- 20) 


Recall that the parameters, (XR, YR, ZR), describe the position of the center element of the 
receive array. 

Inspection of the magnitude plots of these figures clearly show that the computer code 
does indeed handle the different permutations of the X and Y offsets. Consider the 
magnitude plot shown in Fig. 6.10a. Since AY is negative and AX is positive, then from 


the perspective of the receiver looking towards the source, the receiver is located above and 
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Fig. 6.7a Magnitude of the overall system complex frequency response. 
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to the left of the source. Remember also that this geometry relates to the center of the 
receive array. Given this geometry, we expect that if we move left of center, which 
corresponds to moving in a positive direction along the X axis or increasing the element 
index m (see Fig. 6.4), the smaller the acoustic field would become. Conversely, if we 
move to the right of center, the radial distance between the source and the field point 
decreases and the acoustic field increases. Also given this geometry, as we go deeper, 
which corresponds to movement in a positive direction along the Y axis or increasing the 
array element index n, the radial distance between the source and the field point decreases 
and the acoustic field increases. This behavior is clearly demonstrated in Fig. 6.10a. 
Note that in this explanation, it was assumed that the field points did not cross the dividing 
lines where either AX or AY is zero. 

Note the magnitude plot of Fig. 6.7a, and the way the surface caves in by 
approximately 10% around the center element along the Y axis, 1.e.,n = 0. This 
corresponds to the center of the receiver being at the same depth as the source, i.¢., axis- 
axis propagation. This phenomenon is caused by the fact that our expression for the ocean 
medium transfer function is based on the WKB approximation and is not exact. It was 
observed in Chapter 4 that as the propagation vector in the Y direction ky (y) approached 
zero, the approximate WKB solution is not valid since the solution approaches infinity. 
This can also be thought of in terms of the simplifying assumption given by Eq. (4.17) 
being invalid in this region. 

Inspection of Eq. (6.5) reveals that the integration with respect to fr will cause the 
ocean medium transfer function to be calculated for a wide range of values of the 
propagation vector in the Y direction for all spatial locations. Why is it only significant 
when the transmitter and receiver are at close to the same depth? By thinking of the 


problem from the perspective of ray acoustics, with the full-wave solution consisting of an 
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infinite summation of rays, it can be seen that as the transmitter and receiver get closer 
together in depth, the eigenray which connects the two gets closer to being horizontal and 
hence, closer to the invalid region of the WKB approximation. Since the eigenray and, 
those close are a major component of the acoustic signal at the receiver, it is when they get 
close to the invalid region that the effect becomes noticeable. 

This phenomenon also demonstrates an important point about ray acoustics. In ray 
acoustics, rays other than the eigenrays or those in very close proximity are often ignored 
and thought of as being unimportant. As a result, ray acoustics based on the WKB 
approximation would be of no use in an along axis, isospeed propagation problem since the 
resulting acoustic signal would be infinite. However, we have seen that the full-wave 
solution for an invalid region, by considering an infinite set of rays, was still able to give a 
solution to within 10% of the actual value. The full-wave solution is so called because it 
does not make any assumptions about what is important, and integrates over the entire 
region of propagating and evanescent acoustic wave contributions. We have shown that 
this infinite integration allows the solution to be equated to the exact solution for the free- 
space, isospeed problem and gives accurate values for the magnitude and phase of the 


complex acoustic field at a given field point. 


G. VALIDATION OF THE TRANSMIT ARRAY FAR-FIELD DIRECTIVITY 
FUNCTION 

Recall from Chapters 2 and 3 the discussion regarding the far-field directivity function, 
or beam pattern, of a planar array of point sources located in the XY plane. If our 
assumption about the axisymmetric nature of the far-field directivity function is correct, 
then applying the analysis of the preceding sections in this chapter to the overall system 


complex frequency response given by Eq. (3.14) yields 
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1.2(f/co) mr NT 
H(f,r) = 2x f > Ss Cont) Hm tevoey Uncen a | 
0 m=-MT n=-NT 


x erty Vo cf Yn) df,, (6.7) 


which should be a valid approximate full-wave solution for our ocean medium propagation 
problem. 

To validate Eq. (6.7), we will compute the value of the acoustic field at a number of 
field points (with a RLOS = 50 m) due to a 3x3 element transmit array and compare the 
results from exact calculations and the approximate full-wave solution. The transmitted 
electrical signal will be the Lanczos smoothed, rectangular weighted, CW pulse used 
extensively throughout this chapter. Recall that this signal has a maximum frequency 
component of 640 Hz and a carrier frequency of 400 Hz. For the purposes of this 
investigation, we will only consider the carrier frequency of the transmitted signal. 

First we need to check if we are in the far field. The interelement spacing, in both the 
X and Y directions, of the transmit array is equal to one half of the wavelength 
corresponding to the highest frequency component of the transmit signal. Ziomek [Ref. 
3:pp. 118-120] shows that this value for the interelement spacing is a necessary condition 
for avoiding grating lobes under all conditions of beam steering for planar arrays consisting 
of MT x NT (odd) elements. This is the default value used in the computer 
implementation if no user-defined interelement spacing is specified, and is the value for 


which we will conduct this analysis. Combined with Eq. (2.23), this yields the condition 


21 oe 
Ap 


Cc 


RLOS > 





(6.8) 
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which, if satisfied, places the field point in the far field of the transmit array and the far- 
field directivity function is the appropriate description of the beam pattern for this scenario. 
For this case, the right-hand side of Eq. (6.8) is equal to 2.263 m, which places the RLOS 
of 50 m well into the far field. 

If we further assume that the beam pattern is not steered and there is unit amplitude 
rectangular weighting of the array elements, it can be shown that Eqs. (2.28) and (2.30) 


combine to yield 


sin{ur 15/16) sin[vx 15/16] 





PLY) = sacar traGyaneiTl eaiSHGyie? 0) 
where the direction cosines u and v are calculated from 
AX 
= RLOS (6.9a) 
and 
_ AY 
= RLOS - (6.10b) 


Values of the magnitude of the overall system complex frequency response for various 
combinations of AX and AY are detailed in Table 6.6. The exact value of the acoustic field 
is calculated by multiplying the relevant value of the far-field directivity function, given by 
Eq. (6.8), with the exact value of the acoustic field due to a single point source at the same 
RLOS calculated from Eq. (4.25). The simulated value of the overall system complex 


frequency response, on the other hand, is calculated using the approximate-full wave 


12) 


solution given by Eq. (6.7). As can be easily observed from these results, the value of the 
acoustic field calculated using the approximate full-wave solution behaves well and is 
within 0.07% of the exact value for a wide range of direction cosine values. 

Recall from the discussion of the results detailed in Table 6.3 that the numerical 
integration of the overall system complex frequency response was carried out to sufficient 
terms to give a maximum of a 0.01% relative error. Note that the errors detailed in Table 
6.6 are larger than this value. Does this mean that there is an error with the assumption 
that the beam pattern is axisymmetric ? The answerisno. The increased error is a result 
of the manner in which Eq. (6.7) is numerically computed. To aid in the convergence of 
the numerical integration routine, the summations of Eq. (6.7) are taken outside of the 
integral. Therefore, each integration within the summation, with its associated error, 
contributes to the total error of the calculated overall system complex frequency response. 

Next, consider a steered beam pattern. The default value for the direction of 
beamsteering in the computer implementation is the line-of-sight direction between the 
center of the transmit array and the field point. In other words, the beam is steered 
towards the field point. This default value is used if no user-defined direction is specified, 
and this is the case we will investigate for this analysis. Since the beam is steered towards 
the field point, it can be shown that the far-field directivity function given by Eq. (2.28) in 
the direction of the field point reduces to a constant equal to the number of point sources in 
the array. In other words, the array sources coherently add in the direction of the center of 
the main lobe of the beam pattern. As a result, the far-field directivity function of our 
steered transmit array has the value 9 in the direction of the field point. The values of the 
exact and approximate full-wave solutions of the acoustic field due to this steered array are 


detailed in Table 6.7. As can be observed, like the results contained in Table 6.6, the 


2 


exact and approximate results show close agreement, with errors being of the same order of 
magnitude. 

These results confirm that our analysis relating to the far-field directivity function and 
axisymmetric behavior of a planar array of point sources contained in the previous chapters 
was correct. The importance of this result is that it demonstrates that the assumption made 
in some other ocean propagation simulation programs, that axial symmetry can only be 


used for single point sources or vertical line arrays [Ref. 7] is incorrect. 


F. IMPLEMENTATION OF NUMERICAL INTEGRATION 

To implement the integration of Eq. (3.14) with respect to fr, the IMSL10 integration 
routine, DQDAG, was used. This is a double precision numerical integration routine 
which uses a globally adaptive scheme based on Gauss-Konrad rules. Note that the 
integrand of Eq.(3.14) contains both a zero-order Bessel function of the first kind, 
implemented using the double precision IMSL10 function, DBSJO, and complex 
exponential terms. This means that the integrand is oscillatory, becoming more oscillatory 
as frequency and the range between the source and the field point are increased. Asa 
consequence, a Gauss-Konrad rule is used with 30-61 points. 

The oscillatory nature of the integrand means that it is time Consuming to compute 
numerically with a accuracy better than 0.1% - 0.5%. For our computer implementation 
we desired a relative error less than 0.01%, which was not achievable by numerically 
integrating between the appropriate limits, regardless of how much time was allowed. To 
achieve this error criterion, it was necessary to subdivide the original regions of integration 
into 20 subintervals, and to compute each interval separately and sum the result. 

It has been mentioned on numerous occasions that numerical integration is a time 
consuming operation. To illustrate this fact, some typical run times are detailed in Table 


6.8. Note that this program has been implemented in FORTRAN 77 and was run on the 
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AX | Dt (f,fx ) AY Dr (f,fy ) | H(f,r) | | H(f,r) | % error 
( exact ) ( simulated ) 


0 3.0000 Pa) Ze 1.0079 e-2 | 1.0081 e-2 0.020 
20 2.4142 - 20 2.4142 9.2758 e-3 | 9.2795 e-3 0.040 
- 10 2.8478 10 2.8478 1.2907 e-2 | 1.2903 e-2 - 0.031 
-10 2.8478 - 30 1.7654 8.0015 e-3 | 8.0038 e-3 0.029 
35 1.3902 10 2.8478 6.3008 e-3 | 6.3050 e-3 0.067 


Table 6.6 Acoustic field resulting from an 
unsteered, 3 x 3 transmit array. 


AX AY | H(f,r) | | H(f,r) | % error 
( exact ) ( simulated ) 


- 10 10 1.4324 e-2 | 1.4318 e-2 - 0.042 
20 - 20 1.4324 e-2 | 1.4320 e-2 - 0.028 


Table 6.7 Acoustic field resulting from a 3 x 3 transmit array 
steered in the direction of the field point. 
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Naval Postgraduate School IBM 3033/4381 Network. It should be noted that the run 
times shown in Table 6.8 are for a frequency of 400 Hz. Note that the times detailed in 
Table 6.8 are for the simplest case possible. The phase integral term of the ocean medium 
transfer function has a closed-form expression for the isospeed medium case, which means 
that no numerical integration or approximation of the phase integral has to be made. 
Hence, computation time for this component of the integrand is ata minimum. The 
computation time becomes very significant when dealing with receive and / or transmit 


arrays rather than single elements. This is the main reason that most of our analysis has 


been carried out for small values of RLOS. 


RLOS Computation 
(m) Time ( min:sec ) 
10 2:01 
100 2:20 
1000 3:23 
10000 12:18 


Table 6.8 Simulation run times for a single point source 
to a single point receiver problem at 400 Hz. 


125 


VII. COMPUTER SIMULATION RESULTS FOR THE 
SINGLE GRADIENT MEDIUM CASE 


A. EVALUATION OF THE WKB PHASE INTEGRAL 

The results of the previous chapter validated the development carried out in Chapters 2 
through 4 for the case of a free-space, isospeed medium. The next step in the validation 
process is to apply the analysis to a more complex scenario, i.e., the free-space, single 
gradient case. For the purposes of this thesis we will characterize the sound-speed profile 


as being linear, with the speed of sound being given by 
C(y) = Cypert 8(Y - YREF)» (71) 


where g is the sound-speed gradient, and cyprr is the speed of sound at the depth yprr 
The linear approximation of the sound-speed profile was selected as it allows for the easy 
description of a more complex medium by means of linear segments. The aim of this 
chapter is to validate the analysis for a single linear segment before the computer 
implementation is applied to a multiple linear segment sound-speed profile. 

Since the preceding development of the overall system complex frequency response 
given by Eq. (3.14) assumed a free-space problem, all we need to do is apply the condition 


of a single, linear, sound-speed profile to the ocean medium transfer function, given by 


y 
. ne | 
Hite aay) = “yOo) if YOR jomtyWoK-Yn) (4.34) 


2kyYo)|kyoy 
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where 


Zee ne 7 Z 
ky (y)=+2n{(f/ey))"-frl, fr $(f/ey))% (7.2a) 


and 


1/2 
ky (y) =F j2n[ fe -(f/ey))] fe > (£/ely))” (7.2b) 


For the form of the sound-speed profile given by Eq. (7.1), we were unable to find a 
closed form solution to the phase integral in Eq. (4.34). This is in marked contrast to the 
results of the previous chapter, where the phase integral had a closed form expression 
which resulted in the ocean medium transfer function for the isospeed medium requiring no 
numerical integration in its evaluation. 

The method required to evaluate the phase integral of Eq. (4.34) is of critical 
importance. Recall from Chapter 3 that the expression for the overall system complex 
frequency response, given by Eq. (3.14), involved the evaluation of an integral whose 
integrand contained the ocean medium transfer function. In other words, there exists an 
integration within an integration. If a numerical integration routine, such as the IMSL10 
routine, DQDAG, discussed in the previous chapter, is used to evaluate this phase integral, 
then the resulting computational time involved is unmanageable. As an example, the 
computer simulation was run for a single frequency component at a RLOS of 100 m, and a 
single point source and a single point receiver. It was found for this case that the 
maximum allowable processing time of one hour on the Naval Postgraduate School IBM 
3033/4381 Network was attained before a solution could be computed. Consequently, we 


required an alternative method to numerical integration to evaluate the integral. 
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As an approximate solution to the phase integral, we expressed the propagation vector 
components in the Y (depth) direction, given by Eqs. (7.2a) and (7.2b), in terms of their 
binomial series. The reason for doing this was that the integral of each of the terms of the 
binomial series expansion, for both forms of the propagation vector component ky (y), has 
a closed form expression. As a result, the phase integral of the ocean medium transfer 
function can be replaced by a converging series. Taking the integral of the binomial series 


expansion of Eqs. (7.2a) and (7.2b), it can be shown that 
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for the propagating waves described by Eq. (7.2a), and 
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for the evanescent waves described by Eq. (7.2b). The number of terms used in Eqs. 
(7.3) and (7.4) depends on the desired accuracy for the phase integral. Note that this 


method of determining the phase integral is inefficient for small values of ky (y). This is 
due to the fact that for small values of ky (y), f = (f/ c(y)) As a result, the binomial 


expansion, and hence, the phase integral converge very slowly. 
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Recall from Chapter 3 that the plus (minus) sign in Eq. (7.3) is chosen whenever 
Y 2 Yo (y < yo) and corresponds to the field point being deeper (shallower) than the source 
point, 1.e., downward (upward) traveling waves. The minus (plus) sign in Eq. (7.4) 
corresponds to the plus (minus) sign of Eq. (7.3) in order to generate evanescent waves. 
Recall that the full-wave solution can be thought of as an infinite summation of rays. If we 
think of rays from a single omnidirectional point source propagating in a medium 
characterized by a single, linear, sound-speed gradient, it can easily be visualized that some 
rays will never reach the receiver depth while others will reach the receiver depth twice. 
This situation is shown in Fig. 7.1 where it has been assumed that the receiver is below the 
source and the sound-speed gradient is positive. 

This is an important concept when we consider the way in which the overall system 
complex frequency response is calculated. Since the overall system complex frequency 
response is an integral with respect to fr and is a function of frequency, the WKB phase 
integral contained in the integrand is essentially the phase change of a particular ray, 
defined by the values of f and f,, traveling between the source and the receiver depth. For 
each of these rays, we need to decide how it behaves and determine a strategy for 
determining the value of the phase integral. 

This ray path problem was not encountered when dealing with the isospeed medium, 
since for this case the rays travel in straight lines as shown in Fig. 7.2. The situation 
shown in Fig. 7.2 assumes that the receiver is below the source and, as can be seen, all the 
downward traveling waves reach the receiver depth once in an isospeed medium. Asa 
consequence, only the closed form expression of the phase integral for downward traveling 
waves needs to be used in the evaluation of the ocean medium transfer function for 
propagating waves. Similarly, only one form of the phase integral expression needs to be 


used for the evanescent waves. 
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Fig. 7.1 Ray paths in a medium characterized by 
a single, positive, sound-speed gradient. 
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Fig 7.2 Ray paths in an isospeed medium. 
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Inspection of Fig. 7.1 shows the various groups of rays that will be encountered. 
Rays A and B are handled in the same manner. Since neither of these rays is refracted 
back to the receiver depth before reaching the horizontal range of the receiver, the phase 
integral is calculated for the downward traveling portion of the rays. For rays A and B, 
the integral is calculated from the source depth yg meters to the first attainment of the 
receiver depth y meters, i.e., from yg to points 1 and 2, respectively. 

Ray C, unlike rays A and B, is refracted back to the receiver depth before the horizontal 
range of the receiver. As is illustrated in Fig. 7.1, the ray has both downward and upward 
portions that need to be taken into consideration. Consequently, the phase integral is 
calculated for a downward traveling ray from the source to the turning point at point 3, and 
for an upward traveling ray from the turning point back to the receiver depth at point 4. 
This essentially means that the phase integral has been subdivided into two separate 
segments, dependent on which form of the propagation vector given by Eq. (7.2a) is valid. 
Brekhovkikh [Ref. 6:pp. 66-67] also shows that such a ray encounters a the phase shift of 
-m/2 radians at the turning depth, i.e., point 3. After passing through the turning point, the 
upward traveling ray keeps the -7/2 phase shift. 

The final ray type to be considered is that characterized by ray D. As can be seen from 
Fig. 7.1, this propagating ray never reaches the receiver depth since it reaches a turning 
point above the receiver depth at point 5, and becomes an upward traveling ray. To 
evaluate the phase integral, we first consider the downward traveling portion of the ray 
down to the depth of the turning point. We then need to evaluate the integral from this 
point down to the receiver depth. There is clearly no propagating wave component going 
down to the depth of the receiver from the turning point, so any acoustic field reaching the 
depth of the receiver from this point must be evanescent. Accordingly, to evaluate the 


phase integral down to the receiver depth we consider the integrand to be the complex 


131 


propagation vector component given by Eq. (7.4). Recalling the discussion regarding the 
behavior of ray C at a tuming point, this evanescent component can be thought of as what 
gives rise to the -m/2 phase shift in the propagating portion of the ray reflected at the turning 
point. In other words, when the ray goes through a turning point, it gives rise to a phase 
shifted, reflected propagating component and a evanescent component. 

Recall from the previous chapter the importance of evanescent waves. As was done 


for the isospeed medium case, we will use 


1.2 (f/¢9) “a 
Hf,r) = 2n f Hy (f£,,Yo:y) Jo (tft) Ce IO me df. (6.5) 
0 


as our approximate full-wave solution to the ocean propagation problem which means that 
we will be integrating with respect to fr up to a limit of 1.2 (f/co). Since this integration 
limit gave very good results for the isospeed medium case, it is assumed that it will also be 
applicable to the new single, linear gradient problem since the sound-speedprofile is slowly 


changing and does not exhibit any erratic behavior. Note that the region of integration 


(f/c,)<t, <1 2 (tice (7.4) 


describes evanescent waves resulting from the source condition and should not be confused 


with an evanescent wave that results from a ray being refracted through a turning point. 


B. RESULTS 
As discussed in the preceding section, the accuracy of the binomial series 


approximation depends on the number of terms used and for small values of ky(y) the 


number of terms required becomes very large. Just as for the situation in which we tried 
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to evaluate the phase integral using a numerical integration routine, if we specify a 
requirement that the phase integral be accurate to within 0.01% and evaluate the binomial 
expansion approximation for all the values of f used to numerically integrate Eq. (6.5), the 
computation time becomes unmanageable. In this case it was again found that the 
maximum allowable processing time of one hour on the Naval Postgraduate School IBM 
3033/4381 Network was attained before a solution could be computed. 

As a consequence, it was decided to modify the region of integration in the evaluation 
of Eq. (6.5). Instead of integrating over the two regions O<f, $(f/c,))and 
(f/Co) $f, $1.2 (f/Cg), we calculated the approximate full-wave solution of Eq. (6.5) 
using the following modified regions: 0 <f, < LIM (f/c,) for initially propagating waves 
and (1/LIM) (f/cg)<$f, $1.2 (f/c)) for initially evanescent waves. In these 
expressions, LIM is a constant less than unity that basically allows us to trade computation 
time against accuracy. After some investigation of the rate of convergence of the binomial 
series approximation of the phase integral, it was found that solutions could be obtained for 
a scenario in which the RLOS was 100 m, the transmitter was a single point source, and 
the receiver was a 3 x 3 array of point receivers, if the value of LIM did not exceed 0.995. 
It can easily be seen that the closer the value of LIM is to unity, the better our approximate 
full-wave solution should be. Accordingly, we conducted simulations for LIM equal to 
0.99 and 0.995. 

From our observations in the previous chapter, there are two forms of output that are 
useful in determining the performance of our full-wave approximation. First, there is the 
tabulation of the magnitude and phase of the complex acoustic field at the locations of the 
receive array elements. Recall that the usefulness of this form of output was that it allows 
for easy verification that the magnitude of the complex acoustic field is behaving correctly. 


For a given source / receiver geometry, it is an easy matter to determine if the the magnitude 
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of the complex acoustic field is increasing or decreasing according to theory as we go from 
one element to the next. Second, there is the output from the adaptive beamforming 
algorithm, used in the previous chapter, processing the time-domain electrical signal at each 
element of the receive array. This algorithm primarily uses the phase of the received signal 
in its processing. Therefore, it validates how accurate the phase of the predicted complex 
acoustic field is. Consequently, analysis of both these outputs provides a check of the 
accuracy of both the amplitude and phase of the complex acoustic field at a specified field 
point. 

Computer simulation results were generated for two values of the constant LIM, as 
previously discussed, in conjunction with a source / receiver geometry of the form shown 
in Fig. 7.1. The source was placed at a depth of 1,000 m, the receiver was below at a 
depth of 1,050 m, and the RLOS was 100 m._ The sound-speed profile of the ocean 
medium was characterized by the linear sound-speed profile given by Eq. (7.1), where 
YREF= 1,000 m, Cyrgr= 1,475 m/sec and the gradient had the value 0.017 1/sec. The 
aim of investigating these particular values of LIM was to determine whether the small 
region of integration that was discarded was significant. 

The values of the magnitude and phase (in degrees) of the acoustic field for four test 
cases are shown in Tables 7.1 through 7.4. Note that the first value is the magnitude and 
the second value is the phase. Given the placement of the receive array relative to the 
source, i.e., XR = 0 and YR = 1050 m, which is below the source, we would expect the 
magnitude of the acoustic field to decrease with increasing n. This corresponds to reading 
down the table or equivalently increasing the depth of the receive element or field point. 
We would also expect the value of the acoustic field to decrease with increasing absolute 
value of m. This corresponds to reading across the table out from the center or 


equivalently increasing the displacement in the X direction of the element or field point. 
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We expect the magnitude of the acoustic field to decrease in both of these situations since 
their effect is to increase the radial range between the source and the receive element , or 
field point. Note that these are the same results as expected for the isospeed medium case. 
This is because the RLOS is a small value, and hence, the propagating rays have not had 
much of an opportunity to refract. Therefore, a crude approximation is that the medium is 
homogeneous over this short range. 

Case A, shown in Table 7.1, evaluates the overall system complex frequency response 
over the modified region of the initially propagating waves. In other words, we only 
integrate over the region 0 < f, < LIM (f/c,) and ignore the region of initially evanescent 
waves. Case B, on the other hand, shown in Table 7.2, is the same as case A except that 
it also includes the modified region of intitially evanescent waves. Both cases are 
computed for LIM = 0.99. It can be easily observed from these two tables that the 
magnitude of the acoustic field does not behave exactly as expected. In both cases, the 
magnitude of the acoustic field decrease as we read down the table as expected. However, 
in both cases, the second and third rows have the magnitude increasing as we move away 
from m = 0. 

Note that the magnitude and phase of the complex acoustic field is very similar for both 
these cases despite the inclusion of evanescent waves in Case B. At first, this may appear 
to contradict the results of the previous chapter. For the isospeed medium the inclusion of 
evanescent waves had a dramatic effect on the value of the complex acoustic field. So 
what causes this apparently contradictory behavior? The difference between the observed 
results for the isospeed and single sound-speed gradient cases is due to the region of 
evanescent waves considered. The evanescent waves in the region very close to 


f, =(f/C,) are by far the most dominant ones. Accordingly, the modified region of 


evanescent waves ignores the most dominant ones. So why the upper limit of 
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-1 8.17543 e-4 8.18198 e-4 8.17543 e-4 
- 165.452 - 164.725 - 165.452 
0 8.12795 e-4 8.12545 e-4 8.12795 e-4 
139.886 140.616 139.886 
1 7.98127 e-4 7.97160 e-4 7.98127 e-4 
83.940 84.627 83.940 


Table 7.1 Case A: LIM = 0.99, initially evanescent waves not included. 


n/m 1 0 =i} 
-] 8.17524 e-4 8.18179 e-4 8.17524 e-4 
- 165.451 - 164.725 - 165.451 
0 8.12783 e-4 8.12534 e-4 8.12783 e-4 
139.885 140.615 139.885 
1 7.98128 e-4 7.97161 e-4 7.98128 e-4 
83.939 84.626 83.939 


Table 7.2 Case B: LIM = 0.99, initially evanescent waves included. 
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f, =1.2(f/c,)? This upper limit becomes significant when considering very short 
ranges, Le., ranges less than 100 m._ At these short ranges the region of significant 
evanescent waves increases rapidly with decreasing range. The specific upper limit of 
f, = 1.2 (f£/c,) was chosen to give accurate values for the acoustic field calculated from 
the approximate full-wave solution down to a range of 10 m. 

If our reasoning has been correct, then we would expect our values for the computed 
complex acoustic field to improve as we increased the value of LIM closer to unity. Case 
C, which ignores the modified region of evanescent waves, and Case D, which includes 
the modified region of evanescent waves, were both computed for LIM = 0.995 and are 
shown in Tables 7.3 and 7.4, respectively. As can be seen from the values in these tables, 
the same phenomena is demonstrated as was evident when comparing Tables 7.1 and 7.2. 
In fact, it can be easily seen that increasing the region of evanescent waves has had minimal 
effect on the value of the complex acoustic field. Does this mean that our explanaton of 
the importance of the evanescent waves in the region very close to f, =(f/cg) 1s 
incorrect? Let us reconsider the isospeed medium of the previous chapter and see what 
happens when we apply these modified regions of integration. 

Consider the same source / receiver geometry as we have been using for this 
investigation, but apply the approximate full-wave solution to an isospeed problem, where 
the constant speed of sound Co is 1,475 m/sec as was the case in the previous chapter. 
Applying the modified regions of integration, with LIM = 0.999, to the ocean medium 
tansfer function given by Eq. (6.5) results in the values of the complex acoustic field 
shown in Table 7.5. Note the way that the magnitudes of the acoustic field do not behave 
as predicted. In the third row the magnitude of the acoustic field increases as we move 
away from m=(0. Also, note that the magnitudes of the acoustic field increase with 


increasing depth. Both of these results are the opposite from what we expect. Comparing 
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-] 8.29299 e-4 
- 162.942 
0 7.94735 e-4 
142.244 
1 7.61585 e-4 
84.539 


8.28427 e-4 
- 162.348 


7.94784 e-4 
142.804 


7.62560 e-4 
85.126 


8.29299 e-4 
- 162.942 


7.94735 e-4 
142.244 


7.61585 e-4 
84.539 


Table 7.3 Case C: LIM = 0.995, initially evanescent waves not included. 


n/m 1 
-1 8.29110 e-4 
- 162.938 
0 7.94604 e-4 
142.237 
1 7.61599 e-4 
84.529 


8.28241 e-4 
- 162.344 


7.94653 e-4 
142.797 


7.62572 €-4 
85.116 


8.29110 e-4 
- 162.938 


7.94604 e-4 
142.237 


7.61599 e-4 
84.529 


Table 7.4 Case D: LIM = 0.995, initially evanescent waves included. 
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Tables 7.5 and 6.5, we observe that there is a significant difference between the two sets of 
results, yet the only difference between the two was that a modified region of integration 
was used to calculate the values in Table 7.5. Recall that the value of LIM for the values 
calculated in Table 7.5 was 0.999, that is, very close to unity. These results reaffirm our 
earlier conclusions that the evanescent waves in the region very close to f, =(f/Cg) 
contribute very significantly to the overall value of the integral. This implies that around 
this value of fr, there may be a region of relatively constant phase in the integrand of the 
overall system complex frequency response. 

The outputs from the LMS adaptive beamforming algorithm for Cases A through D are 
shown in Table 7.6. Recall the geometry illustrated in Fig. 6.3. The angle definitions are 
sull valid with the exception that they relate to the local angles of arrival of the wave at the 
receive transducer. In other words, for the isospeed medium case the local angle of arrival 
was also the direction of the source since all rays traveled in straight lines. From the 
values for all four cases we see that the range estimates RLOS (X) are in error by an order 
of 10 - 15 % while the estimates RLOS (Y) are in error by up to 65 %. This phenomena 
was also observed for the isospeed case when the effects of evanescent were ignored. For 
that case, 4 was argued that since the overall system complex frequency response is 
axisymmetric, then the relative behavior of the acoustic field in the X direction should be 
minimally affected by the behavior of the acoustic field in the Y direction. Consequently, 
if there is a problem with the description of the full-wave solution in the Y direction, as in 
this case, the value of RLOS (X) should still be fairly accurate, which is the result that we 
observed. 

The angular outputs shown in Table 7.6 agree very closely with theory. Even though 
the angles are the predicted local angle of arrival, for our case where the RLOS is small 


(i.e., 100 m), the effects of refraction will have been small and the local angle of 
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7.33866 e-4 
- 169.589 


7.75856 e-4 
131.862 


8.30442 e-4 
76.103 


7.34792 e-4 
- 168.589 


7.77764 e-4 
132.604 


8.29309 e-4 
76.781 


7.33866 e-4 
- 169.321 


7.75856 e-4 
131.862 


8.30442 e-4 
76.103 


Table 7.5 Case E: LIM = 0.999, initially evanescent included. 


RLOS (X) 
(m) 


est 


RLOS (Y) 


(m ) 


est 


Table 7.6 Output from adaptive LMS beamforming algorithm. 


90.690 
90.688 
115.748 
115.720 


15517 
75.628 
34.343 
34.438 
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arrival should be very close to the angles describing the position of the source relative to the 
receiver. Any effect due to refraction for the given geometry and sound-speed profile will 
be to decrease the angle © slightly and leave ‘Y unaffected. However, our results are too 
inconsistent to tell if there is a noticeable decrease in the value of ©, but they do 
demonstrate that ‘¥ remains unchanged. 

The results obtained from this investigation of a medium characterized by a linear 
sound-speed profile with a single gradient have essentially been inconclusive. It was 
found that the evaluation of the phase integral in the ocean medium transfer function based 
on the WKB approximation is of crucial importance to the performance of the approximate 
full-wave solution. If the linear systems approach, and the coupling equations developed 
in this thesis are to be successfully applied to the WKB approximation of the ocean 
medium, we require a fast, accurate method of determining the value of the phase integral 
to be able to accurately compute the value of the complex acoustic field at some given 
spatial location. 

The method of representing the phase integral in terms of a power series expansion, 
while better that using numerical integration techniques, was inappropriate since there is a 
major contribution to the value of the acoustic field in the region where this expansion 
converges very slowly. The use of the series expansion allowed approximate values of the 
acoustic field to be calculated, whereas the use of the IMSL10 numerical integration routine 


DQDAG took too much computing time and no useful results could be obtained. 
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VII. CONCLUSIONS AND RECOMMENDATIONS 


It appears from our preliminary results that our approach to modelling the ocean 
medium pulse propagation problem, which incorporates both linear systems theory and the 
physics of acoustic propagation, is a valid one. 

A signal generation scheme based on a truncated Fourier series and complex envelope 
representation of narrowband signals was implemented and shown to work correctly. The 
scheme is capable of generating CW and LFM pulses with rectangular, Hamming or 
Hanning amplitude weighting suitable for conversion into a transmitted acoustic signal. 

We have shown that the beam pattern of the transmit array does not violate axial 
symmetry. This is a result of the theory of linear superposition, which allows us in the 
case of our linear system model, to decompose the beam pattern into a sum of 
omnidirectional point sources. 

The model has been extensively tested and validated for the case of a free-space 
isospeed medium. It was found that the formal solution to the propagation problem could 
be modified to include only a small portion of the region of evanescent waves. This 
approximate full-wave solution has demonstrated that it gives extremely accurate results for 
radial ranges between source and receiver down to 10 m._ It was also shown that the 
formal solution was inaccurate in the situation where the source and field point were at or 
very near the same depth. It was concluded that this was due to the nature of the WKB 
approximation rather than a problem with the system model. It is recommended that a 
modified WKB approximation be used in the vicinity of a turning point, which involves 
Airy functions. This modified solution should be investigated in an attempt to determine if 


it improves the model performance for this source / receiver geometry. 
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Inconclusive results were obtained about the model's behavior for an ocean medium 
characterized by a linear sound-speed profile with a single gradient. This was due to the 
difficulties encountered in evaluating the phase integral term of the ocean medium transfer 
function based on the WKB approximation. For this case we were unable to determine a 
closed form expression for the phase integral. Numerical integration techniques were too 
slow and our derived approximation to the integral was unsuited to the numerical behavior 
of the integrand and its effect on the behavior of the approximate full-wave solution to the 
propagation problem. Further investigation needs to be carried out on how to quickly and 
accurately evaluate this integral. Some recommended avenues for further investigation 
ane: 

* to determine whether a different description of the sound-speed profile will yield a 

closed form expression for the phase integral, 


* to find a better approximation of the phase integral than the binomial series 
expansion. 


There is also a question regarding the strategy required to implement the phase integral. 
From a ray acoustics viewpoint, this integral represents the phase change experienced by a 
ray as it travels from the source to the receiver depth. For the inhomogeneous medium, 
these rays will be refracted and there exists the situation were the ray doesn't reach the 
receiver depth at all, or reaches it more than once. How should the integral be evaluated in 
these circumstances? Because of the inconclusive nature of our computer simulation 
results, we were unable to validate or disprove our scheme to compute the phase integral. 

It was demonswated that the outputs from the computer simulation, i.e., the magnitude 
and phase of the complex acoustic field as a function of frequency and spatial coordinates, 
and the time-domain output electrical signal from each element in a receive array were 


implemented correctly. While no matched-field processing was done using this first 
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output, it was demonstrated that the plotting of the magnitude of the complex acoustic field 
was a powerful aid in helping to visualize what the acoustic field looked like. The time- 
domain signals were processed by an external frequency domain adaptive beamforming 
program and the data obtained was very useful in validating the performance of the model. 
This also demonstrated that generated data could be successfully passed to an external file, 
to be processed by an ancillary program of our choosing. Thus, the generated data is 
independent of whatever receive array signal processing we wish to carry out. 

It has been shown that an approximate solution to the linear wave equation can be 
implemented as an ocean medium transfer function. However, except when closed form 
solutions or short recursive computation schemes exist for all terms in the ocean medium 
transfer function, we are constrained by the required computation time. It is recommended 


that the feasibility of the following be investigated to decrease computation time : 


* optimization of the computer code at the expense of program modularity, 
¢ parallel processing techniques, 


* availability of a faster computing system than the IBM 3033/4381 Network. 


The very nature of the full-wave solution to the pulse propagation problem makes our 
model computation intensive. If this method of solution is to be a useful tool in 
investigating realistic scenarios, we need to consider either increasing the speed with which 
the current calculations are handled, or use another technique, such as the Fast-Field- 


Program method of evaluating the wave vector component integrals. 
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